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Abstract 


Multiple  scale  methods,  which  are  based  on  discrete  and  continuous  reproducing 
kernels,  wavelets,  and  integral  window  transforms  are  developed  .  In  this  development,  a 
microscope  is  constructed  with  a  flexible  space-time  localized  window  function  which 
translates  and  dilates  in  space  and  time  to  cover  the  entire  domain  of  interest.  This 
microscope  can  magnify,  examine,  and  record  the  image  of  the  various  scales  and  frequencies 
of  the  response  locally  within  the  support  of  the  window  function.  The  degree  of  magnification 
will  depend  on  the  power  of  the  microscope,  a  flexible  space-scale  and  time-frequency  window 
function.  This  complete  characterization  of  the  unknown  response  is  performed  through  the 
integral  window  transform.  This  localization  process  can  be  achieved  by  dilating  the  flexible 
multiple-scale  window  function.  The  zoom  in  and  zoom  out  capability  of  the  window  function 
is  especially  useful  in  examining  complex  flow  phenomena,  such  as  flow  induced  vibration, 
dynamic  stability  of  flow-structure  interaction,  turbulence  structures,  and  high  frequency 
structural  dynamics  response. 

Introduction 

Traditionally,  researchers  in  computational  analysis  have  concentrated  on  bringing 
more  detail  into  their  structural  system  models,  but  have  ignored  the  multi-time  and  multi- 
spatial  scales  inherent  in  the  structural  response.  To  obtain  resolution  of  medium  frequencies, 
finite  element  meshes  consisting  of  several  thousand  structural  and  continuum  elements  are 
necessary.  To  resolve  the  problem  in  the  time  domain,  either  a  very  small  integration  time 
step  is  required  or  if  modal  superposition  is  used  a  large  number  of  structural  eigenmodes  is 
called  for.  These  solutions  often  require  hours  of  computer  time  on  even  the  latest 
supercomputers.  This  severely  limits  the  usefulness  of  these  computations  since  their  use  in 
the  design  process  is  almost  impossible. 

Conversely,  researchers  in  structural  design  and  analysis  have  concentrated  on  the 
modeling  of  various  structural  components,  usually  via  lumped  parameters,  but  have  often 
used  drastic  simplifications  so  that  the  multi-scales  problem  is  eliminated;  however,  the 
accuracy  of  these  methods  is  limited  to  the  very  low  end  of  the  spectrum. 

In  fluid-structure  stability  problems,  the  unstable  response  often  arises  from  a 
coupling  between  phenomena  associated  with  substantially  different  frequencies.  For 
example,  vortex  formation  in  a  flow  about  a  blade  or  airfoil  is  initiated  by  relatively  high 
frequency  modes  of  the  structure.  The  excitation  forces  generated  by  the  vortex  and  the 
instability  itself  generally  involve  low  modes  of  structural  response.  The  complete  time 
integration  of  the  equations  for  the  fluid  and  structure  combined  can  be  prohibitively  time- 
consuming  in  even  two  dimensions  and  is  beyond  the  capability  of  even  the  largest 
supercomputers  in  three  dimensions.  Thus  it  can  be  seen  that  methods  which  can  effectively 
treat  problems  with  large  ranges  in  scale  are  needed  in  many  types  of  analysis  which  arise  in 
structural  dynamics  problems.  This  report  is  aimed  at  developing  such  methods  and  in 
studying  their  applications  to  a  class  of  structural  dynamics  problems. 
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Abstract 

Multiple  scale  methods,  which  are  based  on  discrete  and  continuous  reproducing 
kernels,  wavelets,  and  integral  window  transforms  are  developed  .  In  this  development,  a 
microscope  is  constructed  with  a  flexible  space-time  localized  window  function  which 
translates  and  dilates  in  space  and  time  to  cover  the  entire  domain  of  interest.  This 
microscope  can  magnify,  examine,  and  record  the  image  of  the  various  scales  and  frequencies 
of  the  response  locally  within  the  support  of  the  window  function.  The  degree  of  magnification 
will  depend  on  the  power  of  the  microscope,  a  flexible  space-scale  and  time-frequency  window 
function.  This  complete  characterization  of  the  unknown  response  is  performed  through  the 
integral  window  transform.  This  localization  process  can  be  achieved  by  dilating  the  flexible 
multiple-scale  window  function.  The  zoom  in  and  zoom  out  capability  of  the  window  function 
is  especiallv  useful  in  examining  complex  flow  phenomena,  such  as  flow  induced  vibration, 
dynamic  staoility  of  flow-structure  interaction,  turbulence  structures,  and  high  frequency 
structural  dynamics  response. 
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1.  Introduction 

In  a  compressible  flow-structure  system,  the  unstable  response  often  arises  from  a 
coupling  between  phenomena  associated  with  substantially  different  frequencies.  For 
example,  vortex  formation  in  a  flow  about  a  blade  or  an  airfoU  is  initiated  by  relatively  high 
frequency  modes  of  the  structure.  The  excitation  forces  generated  by  the  vortex  and  the 
instability  itself  generally  involve  low  modes  of  structural  response.  The  complete  time 
integration  of  the  equations  for  the  compressible  fluid  and  structure  combined  can  be 
prohibitively  time-consuming  in  even  two  dimensions  and  is  beyond  the  capability  of  even  the 
largest  supercomputers  in  three  dimensions.  Thus  it  can  be  seen  that  methods  which  can 
effectively  treat  problems  with  large  ranges  in  scale  are  needed  in  many  types  of  analysis 

which  arise  in  structural  dynamics  problems. 

The  development  of  computational  methods  for  low-frequency  response  3D  structural 
systems  are  now  reasonably  well  established.  On  the  other  hand,  little  success  has  been 
made  in  the  structural  area  characterized  by  multiple  scales,  where  response  is  often 
dominated  by  the  middle  part  of  the  spectrum.  The  author  believes  that  the  multiple-scale 
reproducing  kernel  particle  methods  coupled  with  wavelets  proposed  here  have  great  promise 
in  dealing  effectively  with  the  difficult  structural  response  problems  in  which  the  medium 
frequencies  are  important,  such  as  problems  involving  impact,  dynamic  instability, 
compressible  flow-structure  interaction,  and  other  local  phenomena. 

Because  the  frequency  shift  enables  the  time  step  to  depend  only  on  the  size  of  the 
frequency  band  and  not  on  the  frequency  extent  of  the  load,  these  methods  will  speed  up 
analysis  of  medium  frequency  response  immensely.  This  is  particularly  attractive  for 
structural  response  exhibiting  multiple  time  scales.  This  new  development  will  enable 
engineers  not  only  to  bring  more  detail  into  their  structural  system  models,  but  will  also 
enhance  the  computer  simulations  of  many  classes  of  multi-scale  structural  dynamics 
analysis.  It  will  improve  the  accuracy,  efficiency,  and  reliability  of  dynamic  analysis. 

In  the  next  section,  a  review  of  the  coupled  compressible  flow-structure  interaction  is 
presented.  In  section  3,  the  weak  form  of  the  equations  of  motion  is  given  and  some  currently 
available  solution  methods  are  discussed.  In  section  4a,  the  discrete  orthogonal  reproducing 
kernel  interpolation  functions  are  reviewed.  In  section  4b,  the  multi-resolution  wavelets 
analysis  is  given.  In  section  5,  we  present  our  proposed  approach  to  study  complex 
structures,  the  multiple  scale  reproducing  kernel  particle  methods.  Sample  examples  are 
given  in  section  6. 

2.  Review  of  Coupled  Compressible  Flow-Structure  Interaction 
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The  ability  to  solve  general  classes  of  fluid-structure  interaction  problems  involving 
finite  deformations  and  stability  depends  upon  the  ability  to  solve  the  corresponding 
uncoupled  fluid  and  structural  problems,  and  also  the  ability  to  interface  fluid  and  structural 
subdomains.  During  the  last  two  decades,  considerable  progress  has  been  made  in  the 
solutions  of  free  and  moving  boundary  problems  which  involve  large  fluid  deformations. 
Among  these  methods  are  the  marker-and-cell  (MAC)  methods  (Hirt  [1975,  1983],  Nichols 
and  Hirt  [1971],  Harlow  and  Welch  [1965]),  the  volume  of  fluid  (VOF)  methods  (Amsden 
and  Harlow  [1970],  Harlow  et  al.  [1976]),  moving  mesh  techniques  (Subbiah  et  al  [1989]), 
Eulerian  and  arbitrary  Lagrangian-Eulerian  (ALE)  methods  (Liu  et  al.  [1988,  1991],  Huerta 
and  Liu  [1988],  Belytschko  [1983],  Belytschko  and  Kennedy  [1978]),  the  smoothed  particle 
hydrodynamic  (SPH)  methods  and  the  free  Lagrange  methods  (Monaghan  and  Gingold 
[1983],  Gingold  and  Monaghan  [1982],  Burton  and  Harrison  [1991]).  However,  the  coupling 
of  these  methods  with  deformable  structures  is  not  well  understood  and  often  causes 
difficulties. 

In  the  MAC  method,  a  fixed  or  Eulerian  mesh  is  used  for  the  fluid  calculation  and  a 
Lagrangian  set  of  marker  particles  is  used  to  trace  the  moving  free  surface.  Marker  particles 
can  be  spread  over  all  fluid  occupied  regions  with  each  particle  specified  to  move  with  the  fluid 
velocity  at  its  location.  A  free  surface  is  defined  as  lying  at  the  “boundary”  between  regions 
with  and  without  marker  particles.  Mote  specifically,  a  mesh  is  said  to  contain  a  free  surface 
if  the  mesh  contains  cells  with  markers  and  has  at  least  a  neighboring  cell  with  no  markers. 
The  MAC  method  offers  the  distinct  advantage  of  eliminating  all  logic  problems  associated 
with  intersecting  surfaces.  This  method  is  also  readily  extendable  to  three  dimensional 
computations.  However,  the  MAC  method  requires  significant  increase  in  running  time  and 
storage. 

In  the  VOF  approach,  a  function  F  is  defined  whose  value  is  unity  at  any  point  fully 
occupied  by  fluid  and  zero  otherwise.  Cells  with  F  values  between  zero  and  one  must  then 
contain  a  free  surface.  This  fractional  volume  of  fluid  (VOF)  method  provides  the  same 
coarse  interface  information  available  in  the  marker  particle  method.  In  order  to  trace  the  free 
surface,  a  flow  analysis  network  approach  is  developed.  It  uses  the  flux  calculation  at  the 
boundary  of  each  control  volume  to  check  the  fraction  of  fill  instead  of  tracking  marker 
particles.  One  limitation  of  this  control  volume  technique  is  that  the  time  step  must  be 
controlled  to  ensure  that  free  surface  only  passes  one  control  volume  in  a  time  step. 

In  a  moving  mesh  approach,  a  numerical  grid  generation  scheme,  which  facilitates 
solutions  over  arbitrarily  shaped  boundaries,  has  to  be  developed.  This  approach  numerically 
maps  the  irregular  shape  of  the  flow  field  to  a  more  regular  shape  in  a  computational  domain 
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where  the  governing  equations  are  solved.  The  computation  time  of  this  method  is  relatively 
high  because  of  the  large  number  of  time  steps  and  the  mesh  generation  at  each  time  step. 

There  are  two  approaches  to  ALE  methods.  One  approach  updates  the  solution 
variables  in  a  single  time  step  while  the  other  performs  a  Ugrangian  step  foUows  by  an 
Eulerian  or  remap  (or  advection)  step.  The  latter  strategy  is  often  referred  to  as  an  operator 
split  method.  Similar  to  the  moving  mesh  approach,  a  new  mesh  is  required  for  the  ALE 
calculations.  However,  this  new  mesh  is  required  only  if  the  Lagrangian  mesh  is  too 
distorted.  The  major  advantage  of  the  ALE  methods  is  the  moving  boundaries  can  be 
computed  with  high  accuracy  of  the  Lagrangian  method,  and  the  mesh  can  conserve  its 
regularity  in  avoiding  element  entanglement  Although  a  combination  of  multi-material  ALE 
methods  have  also  been  developed,  the  application  of  this  technique  to  multiple  free  surfaces, 
multiple  materials,  and  fluid-structure  interaction  remains  to  be  explored. 

The  SPH  (smooth  particle  hydrodynamics)  and  the  free  Lagrange  methods  are  based 
on  algorithms  which  are  truly  grid-free  or  mesh-free.  Consequently,  these  techniques  do  not 
have  the  mesh  entanglement  problems  and  at  the  same  time,  maintain  the  high  accuracy  of 
the  Lagrangian  calculations.  Similar  to  a  finite  difference/element  Lagrangian  calculation,  the 
continuous  fluid  is  approximated  by  a  set  of  particles  or  fluid  elements  of  equal  mass.  Unlike 
the  traditional  Lagrangian  flnite  element/difference  calculations,  the  inter-particle  forces 
among  these  smooth  fluid  particles  are  derived  from  the  pressure  and  these  interacting 
pressure  forces  are  governed  by  an  interpolating  function.  Hence,  the  nodal  or  element  forces 
(which  ate  usually  constructed  from  a  flnite  element/difference  mesh)  are  no  longer  needed; 
and  consequently,  these  free  Lagrange  methods  require  no  mesh  and  the  mesh  distortion 
problems  ate  completely  eliminated.  The  energy  equation  is  similarly  defined  for  each  smooth 
fluid  particle.  Since  the  fluid  is  modeled  by  equal  mass  particles,  the  density  for  the  fluid  is 
constructed  by  defining  a  weighting  function  in  which  the  density  of  the  equal  mass  particles 
is  proportional  to  the  number  of  particles  per  unit  volume.  With  these  approximations,  the 
motion  of  the  fluid  is  governed  by  the  movement  of  this  set  of  smooth  fluid  particles,  and  the 
movement  of  the  particles  are  governed  by  the  particle  interacting  forces. 

Although  SPH  methods  work  well  if  there  is  no  boundary  (since  the  boundary  terms 
ate  tossed  out  in  the  formulation,  Libersky  and  Petschek  [1990]),  and  when  the  number  of 
unknowns  (nodes)  is  large;  SPH  methods  are  not  as  accurate  as  the  regular  flnite  element 
methods,  Johnson,  Peterson  and  Stryrk,  [1993].  From  our  study  of  SPH  interpolation  function 
via  a  simple  one  dimensional  (ID)  Galerkin  formulation,  we  found  that  there  is  an  additional 
deficiency  in  the  SPH  formulation.  It  is  related  to  the  boundary  correction  term  of  the 
reproducing  kernel  approximation.  We  shall  make  an  attempt  to  identify  this  deficiency  and 
present  our  view  of  improving  the  SPH  kernel  approximation. 
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After  reviewing  the  moving  least  square  interpolation  functions,  Lancaster  and 
Salkauskas  [1981],  and  the  diffuse  element  methods  (DEM),  Nayroles  et  al.  [1992]. 
Belytschko  et  al.  [1993]  pointed  out  that  an  assumption  made  by  Nayroles  et  al.,  the 
interpolation  coefficients  are  constants,  detracts  from  the  accuracy  of  the  method.  They 
developed  the  Element  Free  Galerkin  Methods  (EFGM)  and  showed  that  by  adding  more 
accurate  derivatives  and  enforcing  boundary  conditions  by  Lagrange  multipliers,  the  methods 
could  achieve  very  high  rates  of  convergence.  From  our  experience,  EFGM  are  more  accurate 
than  the  finite  element  methods,  and  hence,  the  SPH  methods  especially  for  a  small  set  of 
nodes.  One  main  drawback  of  EFGM  is  the  computational  expense,  and  we  found  that  it  is 

more  computationally  intensive  than  the  SPH  methods. 

The  objective  of  the  Reproducing  Kernel  Particle  Methods  developed  by  Liu  et  al. 
[1993],  is  along  the  same  line  of  development  as  the  SPH,  DEM,  and  EFGM  :  to  develop  an 
accurate  and  efficient  mesh  free  interpolation  functions.  A  detailed  discussion  on  smooth 
particle  methods,  the  diffuse  element  methods  and  the  element  free  Galerkin  method,  and  the 
recently  developed  reproducing  kernel  particle  methods,  is  also  given  in  the  paper. 

Since  a  continuous  reproducing  kernel  can  be  derived  for  this  method  and  it  is  also  a 
free  Lagrange  particle  method,  we  shall  label  this  development  as  Reproducing  Kernel 
Particle  Methods  (RKPM).  This  proposed  approach  is  motivated  by  the  theory  of  wavelets 
Chui  [1992]  and  Daubechies  [1992],  in  which  a  function  is  represented  by  a  combination  of 
the  dilation  and  translation  of  a  single  wavelet,  which  is  a  window  function.  In  a  wavelet 
analysis,  similar  to  the  SPH  interpolation  kernel,  the  interpolation  coefficients  are  defined  in 
terms  of  the  integral  window  transform  of  the  window  function  and  the  solution  itself.  In  this 
proposed  study,  we  shall  make  use  of  the  multiple  frequency  bands  and/or  multiple  scales 
properties  of  wavelets  analysis  (multi-resolution  analysis),  and  the  time-frequency  and/or 
space-scale  localization  properties  of  the  continuous  and  discrete  reproducing  kernel 
approximations. 

3.  Weak  Form  of  the  Equation  of  Motion  and  Solution  Methods 

Consider  the  weak  form  of  the  equations  of  motion  which  can  be  written  as: 

K  <5u**,  u*‘>  +  B<8u*‘(xerg),  u*'(xerg)  -  g(x)>  =  f<5uh,  p>  (3.1) 

where  K<-,  •>,  B<-,  •>  and  f<-,  •>  are  the  usual  weak  form  operator  of  the  governing 
equations,  boundary  constraint  operator  on  fg,  and  the  force  assembly  operator,  respectively. 
6u**(x,  t),  u**(x,  t),  p(x,  t)  and  g(x,  t)  are  the  test  function,  trial  function,  body  force,  and 
prescribed  data  on  the  boundary  fg.  The  governing  equations  can  be  those  of  structural 
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dynamics,  fluid  dynamics,  coupled-fluid-structure  interaction  or  structural  acoustics,  and 
among  others.  The  classical  Galerkin  method  is  to  approximate  uh  (which  is  the 
approximation  to  u(x,  t))  by: 

uK*.  0  =  X  C.  iK*  - «.) 


and 


8u‘(x.t)  =  i8C.*(x-«.) 

where  <t)(x  -  xa)  can  be  the  global  finite  element  shape  functions,  Hughes  [1987],  spectral 
functions  Gottlieb  and  Orszag  [1977],  the  smooth  particle  hydrodynamic  (SPH)  interpolation 
kernel  function,  Lucy  [1977]  and  Monaghan  [1988],  multiple  scale  finite  element  functions 
Liu,  Zhang  and  Ramirez  [1991],  or  wavelet-type  bases,  Chui  [1992]  and  Daubechies  [1992], 
etc.  Substitute  Eqs.  (3.2)  into  Eq.  (3.1)  and  solve  for  Ca  for  a  =  1, ...,  A  coefficients  will  result 
in  the  discrete  approximation  of  u(x),  provided  certain  continuity  requirements  are  met 

It  is  noted  that  the  above  solution  procedures  also  hold  for  the  various  approaches  such 
as  the  space-time  discontinuous  finite  elements,  Hughes  and  Hulbert  [1988],  and  Shakib  and 
Hughes  [1991],  deforming  space-time  discontinuous  finite  elements,  Tezduyar,  Behr  and  Liou 
[1992];  the  arbitrary  Lagrangian-Eulerian  (ALE)  and  Eulerian-type  finite  element  methods, 
Liu  et  al.  [1991],  and  among  others.  All  these  methods  employ  the  same  type  of 
interpolations,  Eqs.  (3.2),  except  a  new  set  of  motion,  mesh  motion,  is  introduced  through 
similar  equations  (3.2).  This  additional  motion  is  used  to  control  the  mesh  or  grid 
deformation  so  that  mesh  distortion  can  be  minimized  and  the  nodal  points  connectivity 
through  the  mesh  or  grid  description  would  not  give  negative  Jacobians,  Liu  et  al.  [1988]. 

Equations  (3.2)  give  a  good  approximation  to  u(x,  t)  when  the  measure  of  the  spacing 
among  Xa,  call  it  h,  becomes  smaller  and  smaller.  One  can  then  visualize  the  coefficient  Cg  as 
an  average  value  of  u(x,  t),  and  <j)(x  -  Xg)  is  the  weighting  function.  If  h  approaches  to  zero  Cg 
approaches  to  the  true  solution  at  Xg.  However,  in  practice,  h  does  not  approach  zero  and  if 
u(x,  t)  is  a  very  nonlinear  function,  solving  for  the  average  values,  as  all  of  these  methods  do, 
would  probably  leave  out  the  fine  details  of  the  response  u(x,  t). 

Therefore,  our  objective  is,  instead  of  solving  for  the  averaged  u(xg)  (i.e.,  Cg),  we  shall 
develop  an  alternative  form  of  Eqs.  (3.2),  which  is  based  on  discrete  reproducing  kernels  and 
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integral  window  transforms.  In  this  approach,  if  we  include  the  time  dimension  into  x,  we  can 
think  of  ^(x  -  Xt)  as  a  flexible  space-time  window  function  located  at  Xg.  Since  u(x)  is  an 
unknown  function,  our  goal  then  is  to  interpret  the  coefficient  Cg  as  a  microscope  which 
magnifies,  examines,  and  records  the  image  of  the  response  u(x)  around  Xg.  This  can  be 
achieved  by  employing  a  space-time  localized  window  function  that  can  recover  the  various 
scales  and  frequencies  of  the  response  u(x)  locally  around  Xg.  This  leads  us  to  the  use  of 
scaling  functions  and  wavelets  which  are  discussed  in  the  next  two  sections. 


4a.  Discrete  Orthogonal  Reproducing  Kernel  Interpolation  Functions 

If  <j)(x  -  Xg)  is  chosen  to  be  the  scaling  function,  Chui  [1992],  Cg  is  identified  as  the 
integral  window  transform  of  the  unknown  response  u  and  -  Xg)  over  the  domain.  That  is, 
the  unknown  coefficients  are  to  be  constructed  so  that: 


C,  =  C,(u,  x.)  =  <u,  <!>,>  = 


I 


u(x)  (j>(x-x,)  dx 


(4.1) 


where  <u,  <t>g>  is  the  integral  window  transform,  and  V  is  the  domain  of  interest.  The 
reconstruction  formula  Eq.  (3.2a)  becomes: 


u‘‘(x,  t)  =  X  4>(x*x«) 


(4.2) 


Equation  (4.2)  is  typical  for  a  discrete  reproducing  kernel  Hilbert  space.  In  seeking  for  the 
solution  ub(x,  t)  using  Eq.  (4.2)  and  a  similar  equation  for  5u^(x),  it  is  necessary  to  define 
nodal  point  xj;  nodal  uj  ■  u(xj);  and  particles  with  nodal  mass  AMj  and  nodal  density  pj. 

Consequently,  the  nodal  volume  AVj  is  determined  by  AMj/pj  =  AVj  for  J  =  1 . NP,  where 

NP  is  the  total  number  of  particles  inside  V.  Hence, 

NP  NP 

^  AVj  =  V  ;  and  ^  AMj  =  total  mass  (4.3) 

j=i  j=i 


Using  numerical  integration  in  Eqs.  (4.2)  and  the  above  definitions,  uh(x,  t)  becomes: 

A 

u'>(x,  t)  =  X  <U,  <|>a>  <I>(X-Xa) 

«al 
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rxieB(x,) 

=  X  X  ^(x»-*a)AV,  U,(t) 

where  xj  are  the  quadrature  points,  and  B(x.)  is  the  support  of  Equation  (4.4)  is  a 

discrete  reproducing  kernel  approximation  (DRKA).  since  u^x,  t)  is  interpolated  via  <|)(x-Xa) 
through  the  integral  window  transform  of  u(x,  t).  An  interesting  interpretation  of  Eq.  (4.4)  is 
as  follows.  There  are  a  total  of  A  sampling  windows  spread  over  the  domain  V.  At  each 
sampling  location  Xa.  a  microscope,  which  is  constructed  from  the  integral  window  transform 
<u,  <|)a>,  examines  the  motion  of  particles  AMj  carrying  nodal  values  uj  passing  through  the 
support  B(xa).  If  the  motion  is  very  nonlinear,  the  number  of  "free  Lagrange  particles  AMj 
passing  through  each  sampling  window  <t>(x-Xa)  is  different  from  time  to  time.  Therefore,  Eq. 
(4.4)  maintains  the  attributes  of  the  free  Lagrange  methods,  and  with  the  appropriate  choice 
of  4»(x-xa),  gives  more  accurate  results. 

Upon  examining  Eq.  (4.4),  since  only  a  set  of  nodal  points  or  particles  are  involved  in 
DRKA  methods,  similar  to  SPH  methods,  DRKA  methods  have  no  problems  associated  with 
mesh/grid  distortion  or  negative  Jacobian  of  the  elements  resulting  from  large  deformation. 
Hence,  it  is  suitable  for  large  deformation  and  high  velocity  flow  problems.  However,  unlike 
SPH  methods  in  which  <t)(x-xa)  in  Eq.  (3.2a)  is  the  interpolation  kernel  function  that  mimic  a 
Dirac  Delta  function  and  Ca  is  the  product  of  the  nodal  mass  AMa  and  the  nodal  value  of  the 
response  Ua  (»u(xa))  divided  by  the  nodal  density  pa,  (Monaghan  [1988]).  That  is 

uKx,  t)  =  X  ^  "«  =  X  AV,Ua 

.=i 

Comparing  Eqs.  (4.4)  and  (4.5),  if  only  one  point  integration  is  used  to  integrate  the  integral 
window  transform  <u,  ^a>.  which  is  a  very  bad  approximation  to  the  integral,  especially  for 
very  nonlinear  or  complex  u(x,  t),  the  SPH  and  DRKA  methods  take  a  similar  form.  However, 
if  more  than  one  point  quadratures  are  employed  to  evaluate  the  integral  window  transform 
<u,  <t)a>,  the  difference  between  the  two  methods  is  apparent  It  is  thus  clear  from  Eq.  (4.4) 
why  we  call  Ca  a  microscope  in  which  the  magnification  power  will  depend  on  the  number  of 
particles  examined  under  the  support  B(xa). 

4b.  Multi-Resolution  Wavelets  Analysis 


<t>(x-x,) 


(4.4) 
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In  this  section,  we  wish  to  construct  and  interpret  Ca  as  2i  flexible  power  microscope. 
This  microscope  can  zoom  in  (sharpen  the  window  function)  to  pick  up  the  detailed  structures 
of  the  response  u(x);  and  zoom  out  (widen  the  window  function)  if  no  further  magnification  of 
the  response  is  necessary.  This  zoom  in  and  zoom  out  capability  is  especially  useful  in 
examining  complex  flow  phenomena,  such  as  flow  induced  vibrations;  dynamic  stability  of 
flow-structure  interaction,  turbulence  structures,  structural  acoustics,  high  frequency  structure 
dynamic  response,  strain  localization  problems,  and  other  engineering  disciplines. 

In  this  development,  since  Ca  is  defined  in  terms  of  the  unknown  response  u(x)  and  the 
multiple  scale  window  function  <|>(x-xa),  Eq.  (4.2)  has  to  be  redefined  for  the  multi-resolution 
wavelets  analysis.  The  window  functions  <j>(x-Xa)  will  be  replaced  by  flexible  power  window 
functions  Vma(x);  whereas  the  integral  window  transform  coefficients  Ca  =  <u,  <|>a>,  a  =  1, .... 
A,  become  Cma,  which  are  defined  by  <u,  Tht^  power  of  the  window  function  is 

controlled  by  an  index  m  =  0,  1,  ...,  M.  The  multi-resoluiion  analysis  is  performed  simply  by 
the  dilation  of  the  window  functions.  That  is,  the  discrete  reproducing  kernel  formula  Eq. 
(4.2)  can  be  re-defined  as: 

M  An  M  An 

U''(X)  =  X  S  Vmt>  Vin«(x)  =  X  Z  Cm.Vina(x) 

ib4)  asl  ~  n»>0  (4.6a) 

where  the  dilation  of  Vma(x)  is  defined  through  the  integer  index  m,  and  AO  >  A1  >  ....  >  AM. 
The  arbitrary  scale,  ao  >  0,  is  usually  set  equal  to  2.  Hence,  the  flexible  space-time  and 
scale-frequency  window  functions  are  defined  by: 

¥m.(x)  =  -  ¥(a5“'(x-xj)  a©  >  0  (4  6b) 

and  the  integral  window  transform  <u,  ¥ina>  is  given  by: 

Cm*  =  <u,  ¥iB>^  “  I  ¥in*(x)  u(x)  dx 

A  (4.6c) 

where  V  is  the  space  (or  space-time)  domain  of  interest 

As  can  be  seen  from  Eq.  (4.6b),  the  mother  wavelet  v(x)  can  be  obtained  by  setting  m  = 
0  and  Xa=  0.  Because  of  our  choice  of  indexing,  m  =  0  corresponds  to  the  smallest  window 
width.  The  integral  window  transform  <u,  ¥oa^’  microscope,  can  pick  up  the  very  fine 

scale  and/or  high  frequencies  up  to  the  arbitrary  scale  ao-  For  m  =  M,  it  corresponds  to  the 
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largest  window  width,  and  CMa  =  VMa^*  microscope,  can  pick  up  the  large  scale 
and/or  low  frequency  response. 

It  is  apparent  from  Eqs.  (4.6)  that  the  microscope  Cma  located  at  Xa  can  examine  the 
finest  scales  of  u(x)  simply  by  dUation.  The  magnification  factor  is  determined  by  the  mother 
window  function,  and  the  flexible  space-scale  and  time-frequency  localized  integral  window 
transform  <u,  Vina>- 

Using  a  numerical  quadrature  integration  scheme  to  discretize  Eq.  (4.6c)  gives; 

Cma  =  <U,  <t)ma>  =  X  ¥ina  (*0  AV,  U, 

j  (4.7a) 

where  B™(xa)  is  the  support  of  the  m^  scale  window  function  located  at  Xg  and  AVj  is  the 
nodal  volume  evaluated  at  xj.  Substituting  Eq.  (4.7a)  into  Eq.  (4.6a),  we  obtain  our  desired 
discrete  multi-resolution  wavelet  analysis  of  u(x): 

M  An  |*ieB“(x.)  \ 

u‘'(x)  =  X  X  X  Vma  (Xj)  AVj  Vma  (x)  Uj 

111=0  a=l  V  J  '  (4.8) 

Similar  equation  can  also  be  written  for  5uh(x).  Equation  (4.8)  can  be  written  in  a  more 
familiar  form  of  multi-resolution  analysis: 


M 

u‘‘(x)  =  X  “§>(*)  =  “o(*)  +  ■"  “m(x) 

m=:0  (4.9a) 


It  can  be  seen  from  Eq.  (4.7a)  that  uh(x)  is  a  direct  sum  of  the  M-scale  solution.  For  each 
reflnement  m,  the  m^-scale  interpolation,  which  also  represents  the  m***-frequency/wave 
number  band  of  the  solution,  is  defmed  by: 


An 

US,(X)  =  X  =  Uml(x)  +  U^2(X)  +  •"  +  UiSAm(x) 

a=l 


(4.9b) 


From  Eq.  (4.9b),  we  can  view  each  uS«(x)  as  a  flexible  power  microscope  examining  u(x) 
around  Xa  and  each  uSm(x)  is  interpolated  through  the  window  function  so  that: 
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(4.9c) 


xjeB“(x.) 

Umt(x)  =  X  ' 

i 


and  the  m***  scale-a***  wavelet  shape  function  (NmalU))  of  particle  J  is  given  by  (no  sura  on  ra, 
a  and  J): 

Nm.j(x)  =  ¥m.  (Xj)  AV,  Vn,.  (x)  (4.9d) 


From  Eq.  (4.9c),  the  summation  on  J  for  xj  e  B“(xa)  will  define  the  global  nodal  connectivity. 
However,  unlike  the  usual  finite  e'  ment,  the  sequence  of  numbering  uj  would  not  cause  any 
mesh  distortion  or  negative  Jacobian  problems  since  there  is  no  element 


If  additional  particles  are  added  to  the  domain  (this  can  easily  be  done  by  splitting  up 
particles),  the  adaptive  multi-resolution  analysis  can  be  summarized  as  follows: 

M 

u^Kx)  =s  X  “ni(x)  multi-resolution  analysis 

nt=o  (4.10a) 


^  \ 

=  ^  I  ^  u^  (x)  I  summation  of  all  frequency  bands  of  interest 
in*0  V*»l  / 


(4.10b) 


uj  are  examined  by  the  m*  scale-a* 
wavelet  under  B“(x  J 

xjeB"(x.)  I  M  Am  \ 

=  X  {XX  (Viiu(Xj)AVjVnit(x))  Juj  global  interpolation  functions 
J  InW)  a-l  I 


M  /  An 

=  I  I 

maO  '(si 


Xj£B“(X,) 


X  Nn,^(x)  Uj 


L  J 


(4.10c) 


(4.10d) 


As  can  be  seen,  there  is  no  change  in  the  multi-resolution  analysis  Eqs.  (4.10a)  and  (4.10b), 
since  the  dilation  index  m  and  the  number  of  window  functions  ymafx)  are  the  same.  There  is 
also  no  change  in  the  order  of  the  window  function  as  ¥ina(x)  is  derived  from  ¥oa(x)-  The 
only  change  is  in  the  summation  on  J  under  the  support  B^^fxg)  and  the  resulting 
"nodal/particle”  matrix  wii^  be  larger  because  of  the  additional  global  interpolation  functions  of 
those  nodes/particles.  Similar  procedures  can  also  be  developed  for  the  deletion  of  particles. 
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If  we  can  construct  discrete  multi-resolution  wavelet  functions  according  to  Eqs.  (4.10), 
we  shall  have  a  truly  mesh  or  grid-free  adaptive  method  and  the  interpolation  is  defined 
through  a  set  of  arbitrary-spaced  nodal  points  or  particles.  Moreover,  if  x|/(x)  is  chosen  to  be 
Ck  smooth  functions,  i.e.,  smooth  in  its  derivatives,  the  discrete  multi-resolution 
reproducing  kernel  interpolation  functions  are  also  C^. 


5.  Multiple  Scale  Reproducing  Kernel  Particle  Methods 

The  multi-resolution  wavelet  analysis  given  in  sections  4  is  based  on  an  infinite  domain 
assumption.  To  apply  wavelets  to  analyzing  complex  structures,  this  restrictive  assumption 
is  no  longer  valid.  Another  intrinsic  deficiency  of  orthogonal  wavelet  is  the  stringent 

requirement: 


j  x®v(x)dx  =  0  m  =  0,  l,....q 

where  q  is  the  degree  of  polynomials  of  the  mother  wavelet  v(x).  From  Eq.  (5.1)  it  follows 
that  a  q^^  order  wavelet  can  not  represent  1,  x,  x^, ...,  x**,  parts  of  the  solution.  To  remedy 
these  two  restrictions,  we  shall  represent  a  function  u(x)  by: 

u(x)  =  uw(x)  -i-  P(x)c 

where  u^(x)  is  the  part  of  the  solution  that  is  obtained  by  the  multi-resolution  wavelet 

reconstruction  as  given  in  section  4b.  P(x)=  {Pi(x),  P2(x) . Pn(x)}  and  c  =  {ci,  C2, ...,  Cn}^ 

are  the  vectors  of  the  n  linear  independent  functions  and  unknown  coefficients,  respectively. 
A  superscript  T  denotes  the  transpose.  We  can  consider  the  P(x)c  term  as  the  residual 
representation  of  u(x)  within  a  bounded  domain.  Following  the  procedures  of  deriving  the 
reproducing  kernel  particle  interpolation  functions  (see  Liu  et  al.[1993]),  the  multiple  scale 
reproducing  kernel  interpolation  function  coupled  with  wavelets  can  be  shown  to  be: 


u(x)  =  u^(x)  +  j  C(x,  y,  ao)  a5‘  <|)(^)  u(y)  dy 
-|c(x,  y,  ao)  <t>(^)u'"(y)  dy 


(5.3) 
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where  the  parameter  ao  determines  the  size  of  the  scaling  function  <j)(x).  The  form  of  the 
correction  function  C(x,  y,  ao)  will  depend  on  the  choice  of  P(x).  It  is  noted  that  the  second 
term  in  Eq.(5.3)  is  the  reproducing  kernel  approximation  described  in  Appends  A.  The  first 
term  is  the  multi-resolution  wavelet  part,  whereas,  the  third  term  connects  the  two 
reproducing  methods.  It  is  interesting  to  point  out  that  by  a  proper  choice  of  ao,  the 
contribution  of  the  coupling  term  can  be  shown  to  be  negligible.  With  this  construction,  the 
wavelet  and  the  reproducing  kernel  terms  give  the  high  and  low  frequency  (or  the  fine  and 
coarse  scales)  representations  of  the  solution  u.  It  is  also  noted  that  u^(x)  can  be  expressed 
by  other  continuous  or  discrete  multiple  scale  reproducing  kernels. 

To  further  examine  this  multiple  frequency/wave  number  bands  wavelet  approximation, 
we  let 


M  An 

U'*'(X)  =  X  X  Vm«>  ¥in.(x) 
aW)  •>! 


(5.4) 


Substituting  Eq.(5.4)  into  Eq.  (5.3)  yields: 

M  Mr  if  \ 

u(x)  =  X  X  (  I  ‘‘(y)  Vma(y)  dyU¥nia(x)-¥m.(x)] 

iia=0  as  I  Vv  f 

+j  C(x.  y,  ao)  ai,‘  u(y)  dy 


(5.5) 


where  the  definition  of: 


<u,  Vma>  = 


1 


u(y)  ¥ina(y)  dy 


(5.6a) 


has  been  used  in  Eq.(5.5).  The  approximation  of  the  wavelet  functions  ¥ina(x)  through  the 
reproducing  kernel,  denoted  by  ¥ma(x)  is: 


¥ma(x)  =  C(x,  y,  ao)  a^,*  <|)(^)  ¥ina(y)  dy 


(5.6b) 
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It  is  now  clear  from  Eq.  (5.6b)  that  ?ma(x)  is  simply  an  approximation  of  Vma(x)  via  the 
reproducing  kernel  reconstruction.  It  is  expected  that  Vni.(x)  is  very  close  to  Vma(x)  for  low 
frequency/wave  number  wavelets.  Consequently,  the  contribution  from  the  low 
frequency/wave  number  wavelets  is  close  to  zero.  However,  depending  on  the  choice  of  ao,  <1>. 
and  \|/;  the  reproducing  kernel  might  not  be  able  to  reconstruct  the  high  frequency/wave 
number  part  of  the  solution  (second  term  in  Eq.  (5.5));  and  these  high  frequency/wave  number 
components  can  be  readily  picked  up  by  multi-resolution  wavelets  analysis. 

Presently,  we  are  working  on  the  theoretical  analysis  of  this  type  of  reproducing  kernel 
methods.  Upon  the  understanding  of  Eqs.  (5.6),  we  shall  implement  the  correct  <|),  and  \|/  into 
the  continuous  multiple-scale  frequency  bands  approximation  of  the  response: 

uKx)  =  I  C(x,  y,  ao,  a*")  a^,i  <1>(^)  u(y)  dy  low  frequency  band 

+  X  I  C(x,  y,  ao,  a”)  ao®  v(^)  u(y)  dy  high  frequency  band  (5.7) 

m=l  Jv 

Unlike  orthogonal  wavelets,  equation  (5.7)  holds  for  arbitrary  domains.  Discretization  of  Eq. 
(5.7)  gives  the  desired  multi-grid/multi-resolution  analysis  of  the  complex  dynamic  systems: 

NPO  ,  V 

u‘'(x)  =  X  [C(x,  xj,  ao,  a®)  Axj  uj]  a^^  low  frequency  band 

j=i 

scaling  factor  coefficients 

+  X  [C(x,  xj,  ao,  a”)  Axj  uj]  a-i  v(^)  1st  higher  frequency  band 

j=i 

Ist-scale  wavelet  coefficients 


+ 

NPM  ^  „k  . 

+  X  [C(x,  XJ,  ao,  a”)  Axj  uj]  a-M  M“  higher  frequency  band 

J=i  (5.8) 

Mtb-scale  wavelet  coefficients 
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It  is  noted  that  NPM  ^  Npl  ^  npO;  and  it  constitutes  an  unstructured  multi-grid  analysis. 
We  also  wish  to  emphasize  that  orthogonal  wavelets  are  not  necessary;  hence  there  is  a 
larger  class  of  window  functions  which  can  give  good  time  and  frequency  localization. 

6.  Numerical  Examples 

The  steady-state  advection-diffusion  problem  can  be  stated  for  the  one-dimensional  case 
as  follows: 

u,xx-  au,x  =  b(x)  in 

with  the  boundary  conditions 

u(xg)=ui  onTg 

u,x(xh)  =  u'j  onTh 

where  is  the  domain  (i.e.,  0  ^  x  ^  L  ).  Fg  is  the  boundary  within  essential  boundary 
conditions  and  Fh  the  boundary  with  natural  boundary  conditions.  The  whole  boundary  is 
F  =  FgUFh  and  FgnFh=j0).  (  ),x  denotes  derivatives  with  respect  to  x.  u  is  the  scalar 
unknown,  a  given  constant,  and  b(x)  a  given  source  term.  The  parameter  a  is  the  advetive 
velocity  divided  by  the  diffusion  coefficient  The  Peclet  number  is  therefore 

Pe=a^  (6.2c) 

2 

The  boundary  conditions  are  given.  This  problem  can  be  viewed  as  a  heat  transfer  problem 
with  convective  and  diffusive  heat  transfer.  The  source  term  b(x)  can  be  caused  by  a 
chemical  reaction.  Following  Hughes  et  al.  [9]  the  weak  form  of  equation  (6.1)  with  the  least 
square  term  can  be  written  as 


LI 


(6.1) 


(6.2a) 

(6.2b) 
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(6.3) 


i 


w(u,*x**  u,x-b)  dCl  +  I  (w.xx’S  w,x)  x  (u,xx*^  u,x-b)  dfl  =  0 

JQ 


where  w  is  an  arbitrary  test  function  and  x  is  a  parameter.  Using  the  approximations  u**  and 
w**  for  the  functions  u  and  w,  we  obtain  the  usual  matrix  equation. 


The  following  numerical  examples  use  a  highly  irregular  source  term  to  cause  a 
nonlinear  solution  with  two  peaks.  The  source  term  consists  of  two  terms  that  are  very 
similar  and  summed  together.  Each  part  is 

_  2ciC2tf  1 )  SecKc  t(x-xo))^ 

-  (e2c^x-xo)+i)2 

2c?  Sect(ci(x-xo))^  Tanti(ci(x-xo)) 

(eC3(x-xo)+.er«2(x-xo)) 

c^tf2(x-xo).t.gca(x~-xo))  *  (l.Tanti(ci(x-xo))) 

(eca(x-Xo)+e^j(x-xo))2 

*  (l-TanlT(ci(x-xo))) 

(ecj(x-xo)+e<j(*-*o))3 

I  c2(egi(^-»°>-erc^x-xo))  (l-Tanl<ci(x-xo)))  ^  Ci  Sech(ci(x-xo)f  \ 

(eC3(x-xo)+e^5(x-xo))2  (e«^x-xo^er«^x-xo))  j 


where  the  parameter  xq  governs  the  position  of  the  peak,  ci  controls  the  sharpness  on  the 
right  side  of  the  peak  and  C2  controls  the  decay  on  the  left  side.  In  the  previous  equation  the 
second  indice  for  the  paramters  xq,  ci  and  C2  is  omitted.  Therefore  xq  =  xoi,  ci  =  ch  and 
C2  =  C2i  for  i=l,2. 

The  resulting  source  term  is 


b(x)  =  kibi{x)  +  k2b2{x) 


(6.5) 
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Figure  6. 1  Source  term  used  in  the  examples 

with  the  parameters 


i 

1 

2 

XOi 

5.0 

2.5 

Cli 

50 

50 

C2i 

1.0 

3.0 

_!Si _ 

0.5 

1.0 

where  the  parameter  k  determines  the  size  of  the  peak.  The  source  term  using  these 
coefficients  is  shown  in  figure  6.1. 

The  homogeneous  solution  for  the  advection-diffusion  equation  on  the  domain  0  ^  x  ^  6 
with  boundary  conditions  u(0H),  u(6)  =  1  becomes 

=  (6.6) 

e<^«  -  1 

where  the  particular  solution  with  b(x)  given  in  (6.17) 
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u^x)  =  kiu?|x>fk2u5(x) 


(6.7) 


where 


l-Tanli(ci(x-xo)) 

gC:(x-*o)+erca(*-»o) 


for  i  =  1,2 


(6.8) 


with  the  appropriate  constants  from  the  table  above.  The  parameter  a  is  set  to  a— 1000  in  the 
following  examples. 

The  results  for  the  diffusion  equation  were  obtained  by  setting  the  parameter  a  to  zero 

and  setting  the  boundary  conditions  to  u(0)=0,  u(6)  =  0. 

The  results  for  the  Reproducing  Kernel  Method  are  obtained  by  using  a  Window  function 

W  of  the  form 


W{z)=  e-2'  (6.9) 

and 

ao  =  2jAx'/|^  (6-10) 

To  show  the  influence  of  the  parameter the  solution  of  the  differential  equation  is  calculated 
for  several  j.  Note  that  for  j=-2.2  the  solution  becomes  unstable  and  that  for  a  j<0  the 
solution  approaches  the  finite  element  solution. 

Note  that  the  wavelets  are  scaled  with  a  factor  Ax  to  scale  the  mother  wavelet  to  the  size  of 
the  mesh. 

The  solutions  of  the  Diffusion  Equation  and  the  Advection  Diffusion  Equation  using  several  different 
methods  are  shown  in  Figures  6.1-  6.6  An  error  plot  for  the  Diffusion  Equation  is  shown  in  Figure  6.5.  The 
solution  for  31  nodes  is  not  representative,  the  source  term  is  underintegrated. 
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Figure  6.3  Solution  of  the  Diffusion  Equation 
using  four  Wavelet  dilations 
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31 

61 

121 

241 

Usual  FEM 

0.19545 

0.18977 

0.07468 

0.01827 

0.19511 

0.18954 

0.07829 

0.06641 

RKPMi  =  -1.0 

0.19589 

0.19054 

0.07507 

0.01839 

RKPMi=  1.0 

0.22884 

0.11999 

0.03746 

0.00463 

ERRf!IRSI!l 

0.20792 

0.16593 

0.07138 

IKBBIBBI 

0.19507 

0.19137 

0.07507 

Table  6.1a  Error  of  numerical  solutions  for 
the  Diffusion  Equation 


num.  of  nodes 

31 

61 

121 

241 

Usual  FEM 

0.80144 

0.79926 

0.56435 

0.27424 

RKPM  i  =  -2.2 

0.80144 

0.79910 

0.56383 

0.27401 

RKPMi  =  -1.0 

0.80283 

0.80029 

0.56555 

0.27422 

RKPMi=  1.0 

0.79341 

0.63022 

0.33317 

0.07116 

WiKMBKM 

0.79900 

0.73682 

0.54145 

IBBBSIISIM 

0.80149 

0.79854 

0.55935 

Table  6.1b  Error  in  the  derivative  of  numerical  solutions  for 
the  Diffusion  Equation 
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Figure  6.6  Solution  of  the  Advection-Diffusion  Equation 
using  Gaussian  Reproducing  Kernel  with  j=  -2.2 
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Reproducing  Kernel  Particle  Methods 


Abstract 

A  new  continuous  reproducing  kernel  interpolation  function,  which  explores  the  attractive 
features  of  the  flexible  time-frequency  or  space-wave  number  localization  of  a  window  function,  is 
developed.  This  method  is  motivated  by  the  theory  of  wavelets,  and  it  also  has  all  the  desirable 
attributes  of  the  recently  proposed  smooth  particle  hydrodynamics  (SPH)  methods,  moving  least 
square  methods  (MLSM),  diffuse  element  methods  (DEM)  and  the  element  free  Galerkin  methods 
(EFGM).  The  proposed  method  maintains  the  advantages  of  the  free  Lagrange  or  SPH  methods; 
however,  it  produces  much  more  accurate  results.  It  is  labelled  as  the  reproducing  kernel  particle 
method  (RKPM).  In  computer  implementation,  RKPM  is  shown  to  be  more  efficient  than  the 
DEM,  and  EFGM.  Moreover,  if  the  window  function  is  0“,  the  solution  and  its  derivatives  are 
also  O*  in  the  entire  domain.  Theoretical  analysis  and  numerical  experiments  reveal  the  stability 
conditions  and  the  effect  of  the  dilation  parameter  on  the  unusually  high  convergence  rates  of  the 
proposed  method. 


1.  Introduction 

During  the  last  two  decades,  considerable  effort  has  been  devoted  to  the  development  of 
mesh  free  or  grid  free  interpolation  methods.  In  most  methods,  the  interpolation  functions  are 
usually  established  by  enforcing  certain  continuity  requirements  around  a  set  of  ordered  (equally 
spaced)  points.  However,  due  to  deformation,  this  set  of  points  can  become  highly  disordered  and 
the  accuracy  deteriorates.  In  addition,  if  the  interpolation  methods,  such  as  finite  element  and  finite 
difference  methods,  require  a  mesh  or  a  grid,  the  distorted  mesh  can  terminate  the  calculation  due 
to  mesh  entanglement  problems,  among  others. 
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Among  the  mesh  or  grid  free  methods  are  the  Smooth  Particle  Hydrodynamic  (SPH) 
Methods,  which  are  sometimes  called  the  free  Lagrange  methods.  They  depend  only  on  a  set  of 
disordered  points  or  particles  as  developed  by  Lucy  [1977],  Gingold  and  Monaghan  [1977], 
among  others.  The  Diffuse  Element  Methods  (DEM)  developed  by  Nayroles,  Touzot  and  Villon 
[1992],  and  the  Element  Free  Galerkin  Methods  (EFGM)  recently  proposed  by  Belytschko,  Lu  and 
Gu  [1993]  are  based  on  the  moving  least  square  interpolation  functions  (MLSM)  presented  by 
Lancaster  and  Salkauskas  [1981].  All  these  methods  do  not  require  a  finite  difference  grid  nor  a 
finite  element  mesh.  Furthermore,  if  the  kernel  functions  (used  in  SPH  methods),  the  weighting 
functions  (used  in  MLSM,  DEM,  and  EFGM).  and  their  derivatives  are  continuous,  the  solution 
and  its  derivatives  are  also  continuous.  The  truly  mesh  free,  the  continuous  solution,  and  the 
continuous  derivatives  are  the  key  selling  points  of  these  methods. 

The  most  attractive  feature  of  SPH  methods  in  a  large  deformation  analysis  is  the  free 
Lagrange  concept.  Although  SPH  methods  work  well  if  there  is  no  boundary  (since  the  boundary 
terms  are  tossed  out  in  the  formulation,  Libersky  and  Petschek  [1990])  and  when  the  number  of 
unknowns  (nodes)  is  large,  SPH  methods  are  not  as  accurate  as  the  regular  finite  element  methods, 
Johnson,  Peterson  and  Stryrk,  [1993].  From  our  study  of  the  SPH  interpolation  function  via  a 
simple  one  dimensional  (ID)  Galerkin  formulation  (see  the  numerical  examples.  Section  7),  we 
found  that  there  is  an  additional  deficiency  in  the  SPH  formulation.  It  is  related  to  the  boundary 
correction  term  of  the  reproducing  kernel  approximation.  We  shall  make  an  attempt  to  identify  this 
deficiency  and  present  our  view  for  improving  the  SPH  kernel  approximation. 


After  reviewing  the  moving  least  square  interpolation  functions  and  the  diffuse  element 
methods,  Belytschko  et  al.  [1993]  pointed  out  that  an  assumption  made  by  Nayroles  et  al.  [1992], 
the  interpolation  coefficients  are  constants,  detracts  from  the  accuracy  of  the  method.  They  showed 
that  by  adding  more  accurate  derivatives  and  enforcing  boundary  conditions  by  Lagrange 
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multipliers,  the  method  could  achieve  very  high  rates  of  convergence.  From  our  experience, 
EFGM  are  more  accurate  than  the  finite  element  methods,  and  hence,  the  SPH  methods  especially 
for  a  small  set  of  nodes.  One  main  drawback  of  EFGM  is  the  computational  expense,  and  we 
found  that  it  is  more  computationally  intensive  than  the  SPH  methods. 

The  objective  of  this  paper  is  along  the  same  line  of  development  as  the  SPH,  DEM,  and 
EFGM  :  to  develop  accurate  and  efficient  mesh  free  interpolation  functions.  Since  a  continuous 
reproducing  kernel  can  be  derived  for  this  proposed  method  and  it  is  a  free  Lagrange  particle 
method,  we  shall  label  this  development  as  Reproducing  Kernel  Particle  Methods  (RKPM).  This 
proposed  approach  is  motivated  by  the  theory  of  wavelets  (Chui  [1992])  where  a  function  is 
represented  by  a  combination  of  the  dilation  and  translation  of  a  single  wavelet,  which  is  a  window 
function.  In  a  wavelet  analysis,  similar  to  the  SPH  interpolation  kernel,  the  interpolation 
coefficients  are  defined  in  terms  of  the  integral  window  transform  of  the  window  function  and  the 
solution.  We  shall  borrow  three  key  ideas  from  wavelet  analysis:  the  integral  window  transform, 
the  dilation  and  translation  of  a  window  function,  and  the  continuous  and  discrete  reproducing 
kernel  i^)proximations.  It  is  noted  that  the  window  functions  used  in  this  paper  are  not  wavelets. 
Good  candidates  for  the  window  functions  are  the  scaling  functions  used  to  produce  wavelets  since 
the  scaling  functions  can  be  constructed  to  be  orthogonal  with  respect  to  its  translates. 

We  shall  show  the  similarities  between  the  smooth  particle  hydrodynamic  methods,  the 
diffuse  element  methods,  the  element  free  Galerkin  methods  and  the  reproducing  kernel  particle 
methods.  We  shall  also  show  that  SPH  and  RKPM  are  indeed  developed  through  a  continuous 
reproducing  kernel  approximation;  whereas,  DEM  and  EFGM,  like  finite  elements,  are  developed 
through  a  discrete  reproducing  kernel  approximation.  As  a  by  product  of  this  development,  the 
concept  of  the  dilation  of  a  window  function  will  be  used  to  explain  why  the  accuracy  of  the 
diffuse  element  methods  decreases  relative  to  the  element  free  Galerkin  methods. 
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In  the  next  section,  some  preliminary  concepts  of  integral  window  transform  and  SPH 
interpolation  kernel  functions  are  reviewed.  In  Section  3,  the  reproducing  kernel  particle 
interpolation  functions  are  derived.  In  Section  4,  the  effect  of  the  dilation  parameter  on  the 
reproducing  kernel,  time-frequency  or  space-wave  number  locaUzation,  and  the  stability  condition 
is  discussed.  In  Section  5,  some  examples  of  the  reproducing  kernel  window  functions  are 
presented.  The  similarities  among  SPH,  DEM,  EFGM.  and  RKPM  interpolation  functions  are 
given  in  Section  6.  Numerical  experiments,  which  confirm  the  theoretical  analysis,  are  presented  m 
Section  7,  followed  by  a  conclusion. 


2.  Preliminaries 

2.1  Dilation  and  translation  of  a  window  function 


Let  X  denote  the  spatial  coordinates.  If  Ofx)  is  a  window  function  located  at  x  =  0,  which 
has  a  support  of  6(x),  then 


<I)(x)  0  in  B(x) 


(2.1a) 


d)(x)  =  0  outside  B(x)  (2.  lb) 

The  dilation  and  translation  of  d>(x),  denoted  by  d>ab(x),  is  defined  as: 

«I>ab(x)  =  E(a)  <D(^^)  a>0  (2.2) 

is  a  window  function,  located  at  x  =  b  with  a  support  scaled  by  the  dilation  parameter  a.  The 
function  E(a)  appearing  in  Eq.  (2.2)  scales  Oa5(x)  such  that: 
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(2.3) 


<I>ab(x)  dR,  =  <t>(x)dR,=  l 

/R,  Jut 

when  the  support  B(x)  is  within  the  spatial  region  of  interest,  Rx-  It  is  noted  that  when  b  is  close 
to  the  boundary  of  Rx.  9Rx.  die  integral  of  I^x  '^dl  be  less  than  1.  We  believe  that 

this  is  a  drawback  for  the  SPH  methods,  as  well  as  all  other  reproducing  kernel  methods,  including 
wavelets,  which  assume  the  region  is  unbounded. 

2.2  Integral  window  transform  and  SPH  interpolation  kernel  functions 

The  integral  window  transform  of  a  real  function  u(x)  with  a  real  window  function  <bab(x) 
is  defined  as: 


<  U  ,  <I)ab  )  = 


E(a)  <I)(^)u(x)  dR, 


(2.4) 


As  a  matter  of  fact,  one  of  the  main  concepts  of  the  SPH  method  is  to  find  a  suitable  smooth 
reproducing  kernel  function  <I>(x)  that  mimics  the  Dirac  Delta  function.  Hence,  when  a  is  chosen  to 
approach  zero  and  when  b  =  x,  the  reproducing  kernel  approximation  of  u(x),  denoted  by  u^(x), 
is  given  by: 


uHx)  =  <  u  ,  <I>ax  >  = 


E(a)<D(^)u(y)  dRy 


(2.5) 
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Discretizing  the  integral  of  Eq.(2.5)  by  NP  distinct  nodes  (points)  using  a  numerical  quadrature 
formula  gives: 


uKx)  s  £  E(a)  u(xj)  AVj  (2.6a) 

j=i 

Equation  (2.6a)  can  be  written  in  a  more  familiar  notation;  in  terms  of  generalized  global  shape 
functions  Nj(x): 


uKx)  =  £  Nj(x)  uj  ;  Nj(x)  =  E(a)  AVj  no  sum  on  J  (2.6b) 

J=i 

where  uj  s  u(xj)  and  A  Vj  0  is  the  nodal  domain  (volume  in  three  dimensions  (3D),  area  in 
2D,  and  length  in  ID)  associated  with  quadrature  points  xj.  The  sum  of  all  AVj  gives  the  total 

dtxnain  V.  That  is: 


NP 

X  AVj  =  V  (2.6c) 

j=i 

The  SPH  methods  use  a  similar  interpolation  formula  as  given  in  Eqs.  (2.6).  Instead  of  using  AVj 
and  E(a)  4>(*f^--),  AMj  and  pj  E(a)  d)(^^^  -)  are  employed.  AMj  and  p;  are  the  Jth  particle 

mass  and  density,  respectively.  With  this  substitution  in  Eq.(2.6a)  or  (2.6b),  the  SPH 
approximation  can  be  used  in  standard  interpolation  Galerkin,  collocation,  or  spectral  methods,  but 
the  particle  methods  use  information  from  a  set  of  disordered  points  based  on  kernel  estimation, 
Monaghan  [1982].  It  is  also  pointed  out  by  Monaghan  [1988],  SPH  works  quite  well  for 
arbitrarily  moving  fluid  if  the  number  of  particles  is  large  and  in  the  absence  of  boundaries; 
however,  there  is  not  a  systematic  way  to  handle  moving  fluid  with  rigid  or  moving  boundaries  in 
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SPH.  Furthermore,  it  is  not  clear  how  to  generalize  SPH  to  non-uniform  mass  particles  or  the 
effect  of  the  dilation  parameter  "a"  on  the  accuracy  of  the  solution. 

2.3  Moments 

We  shall  define  the  following  moments  for  the  window  function  d>(x). 

mo(x)  =  I  p  <I>(y)  dRy  zero  moment  (2.7a) 

mi(x)  =1  p  Yi  d>(y)  dRy  first  moment  (2.7b) 

Jkm) 


mij(x) 


P  yi  Ti  <I>(y)  dRy 


cross  moment 
i=j  second  moment 


(2.7c) 


In  the  above  equations,  the  integral  is  evaluated  with  respect  to  the  support  B(x).  Hence,  if  B(x)  is 
close  to  the  boundary  of  the  spatial  region  R^,  mglx)  is  less  than  one.  The  subscripted  indices  i 

and  j  take  values  from  1  to  NSD,  where  NSD  is  the  number  of  space  dimensions.  If  ^(y)  is 
symmetric,  mi(x)  =  0  in  the  interior  of  Rx  and  mjfx)  ^  0  when  x  is  close  to  the  boundary,  rnjjfx) 
(no  sum  on  i)  doiotes  the  second  moment  of  <I)(x)  in  the  \[  direction;  and  mij(x),  i  ^  j,  denotes  the 
cross  moment  Definitions  for  higher  order  moments  can  be  defined  in  a  similar  fashion.  We  shall 
employ  these  moments  to  analyze  the  reproducing  kernel  particle  interpolation  functions  which  is 
described  next 
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3.  Reproducing  Kernel  Particle  Interpolation  Functions 


The  objective  is  to  use  the  concept  of  reproducing  kernels  and  the  local  character  of  the 
window  function  to  develop  an  accurate  reproducing  kernel  function  from  a  suitable  smooth 
window  function  <I>ab(x)  multipled  by  a  correction  function  C(a,  x,  b).  If  both  <I>ab(x)  and 
C(a,  X,  b)  are  smooth  functions  within  the  spatial  region  Rx,  that  is,  the  functions  and  its 
derivatives  are  continuous,  then  we  have  developed  global  interpolation  functions  that  do  not 
require  a  finite  element  mesh  nor  a  finite  difference  grid.  In  particular,  unlike  the  SPH  methods,  the 
dilation  parameter  a  can  take  a  fairly  large  range  of  positive  values  provided  certain  stability 

conditions  are  met. 

Our  goal  then  is  to  develop  reproducing  kernel  particle  interpolation  functions  which  will 
have  the  following  merits. 

1)  If  d>ab(x)  is  an  even  function,  the  correction  function  C(a,  x,  b)  should  be  one  when  the 
support  of  <I>ab(*)  is  not  close  to  the  boundary.  8Rx;  it  differs  from  one  when  <I)ab(x)  is  dose  to 

the  boundary. 

2)  A  truly  element  or  mesh  free  particle  method  similar  to  SPH  methods  with  much  better 
accuracy,  especially  when  the  number  of  particles  is  stmll. 

3)  Similar  to  MLSM,  DEM,  and  EFGM,  RKPM  provides  smoother  approximations  of  the 
solution  as  well  as  its  derivatives;  however,  RKPM  is  computationally  more  efficient  and  a 
mathematical  analysis  of  the  RKPM  is  also  available. 
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Consider  a  function  u(x)  represented  in  terms  of  m  linearly  independent  functions  Pi(x) 


such  that: 


m 

U(X)  =  X  Pi(x)  Ci 
i=l 


(3.1a) 


or  in  matrix  notation: 


u(x)  =  P(x)  c  (3.1b) 

where  P(x)  =  {Pi(x),  P2(x),  ...  .  Pm(x)},  a  vector  of  m  linear  independent  functions,  and 
c  =  (ci,  C2, ... .  Cm}^  are  the  unknown  coefficients.  A  superposed  "T"  denotes  the  transpose.  In 
order  to  define  c  in  terms  of  the  solution  locally  around  any  point  x,  we  multiply  both  sides  of 
Eq.(3.1b)  by  p(y)  P’^(y)  and  perform  the  integral  window  transform  with  respect  to  a  positive 
even  window  function  <I>ax(y)  yield  (Note:  x  has  been  replaced  by  y  in  Eq.(3.1b)): 

<  pP’ru  ,  ^>ax  )  =  (  pP’^P  .  <l>ax  >  c  (3.2a) 

or  the  vector  of  coefficients  c  is  solved  in  terms  of  the  solution  u: 

c  =  M‘‘(x)  (  pP^u  ,  <l>ax  >  (3.2b) 

where  the  m  X  m  non-singular  matrix  M(x)  is  denoted  by: 

M(x)  =  (  ppT'P  .  (Dax  >  =  I  p(y)  P’^(y)  P(y)  E(a)  <D(^)  dRy  (3.3) 
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It  is  noted  that  M(x)  is  a  continuous  function  of  the  translation  (x)  of  the  window  function  d*ax(y)* 
Substitute  Eq.(3.2b)  into  Eq.(3.1b)  gives  the  approximation  of  u(x),  denoted  by  uh(x)  through  a 
continuous  reproducing  kernel: 

uHx)  =  (  P  P(x)  M-I(x)  U  ,  Oax  >  (3-4) 

Using  the  definition  of  the  integral  window  transform,  u^fx)  can  be  shown  to  be: 

u^fx)  =  I  k(a,  X,  y)  u(y)  dRy  (3.5) 

Jrx 

where  the  reproducing  kernel,  which  is  a  modified  window  function,  is  : 

k(a,  X,  y)  =  kax(y)  =  C(a,  x,  y)  E(a)  (3.6a) 

and  the  function  C(a,  x,  y)  is  given  by: 

C(a,  X,  y)  =  p(y)  (  P(x)  M-‘(x)  P'^(y)  )  (3.6b) 

To  write  Eq.(3.5)  in  a  discrete  reproducing  kernel  particle  form,  the  integral  of  Eq.(3.5)  is 
discretized  by  NP  distinct  points  using  a  numerical  quadrature  formula  to  yield  the  usual 
approximation  formula: 


NP 

u''(x)  =  X  Nj(x)  Uj 


J=i 


(3.7a) 
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and  the  reproducing  kernel  particle  interpolation  functions  are  given  by; 


Nj(x)  =  C(a.  X.  xj)  E(a)  AVj  (3.7b) 

Comparing  Eq.(3.7b)  with  the  SPH  interpolation  formula.  Eq.(2.6b),  we  believe  that  the 
discretized  correction  function 

C(a,  X,  xj)  =  p(xj)  (  P(x)  M-‘(x)  pTCxj)  )  (3.7c) 

will  improve  the  accuracy  of  the  interpolation  kernel  tremendously,  partly  due  to  boundary 
corrections.  Depending  on  the  choice  of  P(x),  we  shall  show  that  the  correction  function  is 
composed  of  the  moments  defined  in  the  previous  session.  In  particular,  if  P(x)  is  chosen  to  be 
constant  and  linear  polynominals: 


P(x)  =  { 1,  xi,  X2,  X3}  (3.8a) 

it  is  shown  that  the  continuous  function  C(a,  x,  y)  takes  the  following  form: 

C(a,  X.  y)  =  Ci(a.  x)  +  C2(a.  x).(^^-^)  (3.8b) 

In  Eq.(3.8b),  Cj  and  C2  are  continuous  scalar  and  vector  functions  of  the  zeroth,  first  and  cross 
moments.  A  dot  denotes  the  inner  product.  It  will  also  be  shown  that  Ci(a,  x)  =  1  and 
C2(a,  x)  =  0  when  the  support  of  ‘I’ax^y)’  not  close  to  the  boundary  of  the  domain; 

whereas  Ci(a,  x)  ^  1  and  C2(a,  x)  ^  0  when  the  window  function  moves  close  to  the  boundary  of 
the  domain.  It  is  noted  that  M(x)  needs  to  be  computed  only  once,  even  in  a  large  deformation  free 
Lagrangian  or  Eulerian  analyses.  In  the  latter  case,  the  density  p  is  chosen  equal  to  one. 


II 


Remark  :  It  can  be  seen  that  if  u(x)  =  P(x),  then: 


i' 


P(X)  =  I  p(y)  P(x)  M-I(x)  p%)  E(a)  P(y)  dRy 

=  P(x)  M'Ux)  I 

=  P(x) 


p(y)  P’^(y)  P(y)  E(a)  dRy 


(3.9) 


because  of  the  definition  of  M(x).  Therefore  if  P(x)  is  chosen  to  be  { 1,  xi.  X2.  X3.  xi^, the 
reproducing  kernel  will  satisfy  the  usual  isoparametric  finite  element  properties.  In  a  discrete 
approximation,  we  have: 

NP  NP 

^  Nj(x)"=  1  ; 

j=i 


4.  Effect  of  the  Dilation  Parameter  on  the  Reproducing  Kernel  and  Stability 
Condition 

In  this  section,  we  restrict  our  discussion  to  one  dimension.  To  explore  the  accuracy  of  the 
RKPM  interpolation  functions,  we  shall  relate  the  dilation  parameter  a  (sometimes  called  a  scale)  to 
the  frequency  content  (in  time);  or  the  wave  number/wavelength  content  (in  space).  For  simplicity, 
we  let  X  be  the  time  axis  and  (O  is  the  frequency  axis.  A  similar  interpretation  can  be  made  when  x 
and  (0  are  the  space  and  the  wave  number,  respectively. 

4.1  Frequency  band 
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Following  Chui  [1992],  we  define  the  center  x*  and  radius  Ai])  of  a  window  function  <I>(x) 


by: 


X*  = 


‘i>(x)F 


dx 


(4.1a) 


A(j>  = 


l*IL 


1/2 


(X-X*)‘ 


|<I>(x)P 


dx 


(4.1b) 


The  width  of  the  window  function  <t>(x)  is  defined  by  2A(j).  The  norm  of  ^>(x)  is  defined  as: 


|(d||2  =  <<I),(I))‘'^ 


(4.2a) 


For  our  choice  of  the  scaling  parameter,  the  norm  of  ^>ab(x)  is  related  to  ‘I>(x)  by: 


(4.2b) 


Suppose  that  <I>(x)  is  any  function  such  that  both  <l>(x)  and  its  Fourier  transform  <I>(CD)  are 
window  functions.  Following  the  definitions  Eqs.(4.1),  the  center  frequency  O)*,  and  the  radius  of 

A<|>  of  4K(i))  we  given  by: 


0)* 


(0 1  <I>((0) 


do) 


(4.3a) 
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(4.3b) 


4.2  Time-freqency  band  of  the  parent  reproducing  kernel  function 

The  center  and  radius  in  Section  4. 1  are  defined  over  an  unbounded  domain  so  that  x*  and 
A(j)  are  invariant  with  respect  to  the  translation  (x).  Since  the  kernel  function  is  a  function  of  y  and 
it  is  defined  within  R*  only,  we  redefine  a  parent  kernel  function  as: 

k(y)  =  C(y)  a-i  d>(y)  (4.4) 

where  C(y)  is  equivalent  to  C(a.  x,  y)  of  Eq.(3.6b)  and  the  arguments  a  and  x  are  dropped  here 
because  they  are  implicit  parameters  and  only  y  is  a  variable.  Let  us  denote  the  center  and  radius  of 
k(y)  as  X  and  Ak,  respectively. 


■  =  — L  j  y|c(y)a-i 

llkigA, 


,-i  cb(y)Pdy 


(4.5a) 


Ak  = _ i — 

liklb 


I 


(y-x)Hc(y)a-i  cl)(y)|^dy 


1/2 


(4.5b) 


If  the  Fourier  transform  of  k(y)  is  denoted  by  k(co),  the  center  (co)  and  the  radius  (Ak)  of 
k(Q))  in  the  frequency  space  are  given  by: 
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CO 


Icill  I 


co[k(co)|  dco 


(4.6a) 


Ak  =  I  J  (CO  -  co)^  I  k(co)  I  ^  (ico 


(4.6b) 


From  the  parent  kernel  function  k(y).  the  two- parameter  reproducing  kernel  fucntion 
kax(y)  =  k(a,  x,  y)  can  be  generated  by  translation  (x)  and  dilation  (a): 


lc..(y)  =  C(3^)a  ' 


(4.7) 


The  relationship  between  the  norms  of  k(y)  and  kax(y)  arc  obtained  by: 


l|k«lb=|[  lc(4^)a-'<t(i^)Pdy 


1/2 


=  a‘'2 


j  |c(z)a-‘<I>(z)p 


dz 


1/2 


(4.8a) 


(4.8b) 


=  ai/211k|b 


(4.8c) 


4.3  Time>frequency  or  space-wave  number  localization 

Because  of  the  linear  translation  of  the  reproducing  kernel  k^^fy).  the  integral  window 
transform  of  the  response  u  with  kg^^fy): 
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<U,  kax>  =  C(^)  a-i  <D(^)  u(y)  dy 

Jf^x 


(4.9a) 


localizes  the  response  with  a  time  window  or  a  space  window: 


[  X  +  ax  -  aAk  ;  x  +  ax  +  aAk  ] 


(4.9b) 


Furthermore,  since  the  domain  is  bounded  ,  x*  and  Ak  have  been  replaced  by  x  and  Ak  such  that: 


X  =  X*  =  constant  ;  Ak  =  Ak  =  constant  ifB(x)in  R*  (4.9c) 


X  =  X*  =  variable  ;  Ak  =  Ak  =  variable  if  B(x)  is  close  to  Rx  (4.9d) 


In  Eq.(4.9d),  x*  and  Ak  have  to  be  computed  according  to  Eqs.(4.5).  The  center  and  radius  of  the 
reproducing  kernel  function  k„(y)  are  shown  to  be: 


radius: 


- - Lf  y|c(i^,a-.4.(i^)Pdy 

=  — 1 — I  (x  +  az)|  C(z)a’' iWz)  pdz  =x  +  ax 

llklgj,. 

[y-(x  +  ax)]2lC(i^)a-i*(i^)fdy 


(4.10a) 


1/2 


I 


= .  .  -  n  (z  -  x)^  I  C(z)  a-i  <I>(z)  I'^dz  )  =  a  Ak 

k  1 12 


(4.10b) 
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The  window  width  is  defined  as  2aAk.  Therefore,  the  integral  reproducing  kernel  window 
transform  can  be  interpreted  as  a  time  or  space-bcalization  of  the  response  (Chui  [1992]).  Hence, 
the  integral  window  transform  Eq.(4.9a)  gives  local  information  of  u  with  the  time  or  space 
window  in  Eq.(4.9b). 

Since  the  Fourier  transform  of  the  reproducing  kernel  k^j^Cy)  denoted  by  kax(ta)  is  also  a 
window  function,  kax(co)  is  given  by: 

kax(®)  =  I  e-‘“y  {C(5^)  a-i  dy 


g-i(aoi)z  {c(z)  a**  <I>(z))  dz  =  a  e’'“*  k(a  to) 


(4.11) 


To  study  u  in  the  frequency  domain,  we  shall  employ  the  Parseval  identity  (Daubechies[1992])  in 
which  the  inner  product  of  any  two  continuous  functions  is  related  to  the  inner  product  of  the 
Fourier  transform,  denoted  by  u  and  g  through; 

<u.  g>=^(u.  g>  (4.12) 


Therefore,  if  we  let  g(y)  =  k„(y)  ,  the  integral  reproducing  kernel  transform  Eq.(3.5)  becomes 


Equation  (4.13)  reveals  that  the  integral  reproducing  kernel  window  transform  also  gives  local 
information  of  u((0)  with  sl  frequency-window 


ja  .  •  fa  +  Ak 

.a  a  ’  a  a  J 


(4.14) 


and  a  bandwidth  equal  to  2Ak/a  .  It  is  noted  that  the  ratio  of  the  center  frequency  and  the 
bandwidth  is  equal  to  m/(2Ak),  which  is  independent  of  the  scaling  parameter  when  the  support  of 
the  window  function  is  within  the  domain  Rx  We  wish  to  construct  C(y)  =  C(a,  x,  y)  so  that  the 

ratio  is  fairly  constant  when  x  is  closed  to  the  boundary.  Employing  the  definition  of  the  center 
frequency,  Eq.(4.6a),  and  the  linear  transform  (o'  =  a(o,  the  center  frequency  of  k(a  (o)  is  shown 

to  be: 


center  =  ■ 


k(ao)) 


£ 


(0  k(a  CO)  dco 


(O’ I  k((o’)  p  d(o’  = 


(4.15a) 


Similarly,  using  Eq.(4.6b),  and  the  linear  transformation  co'  =  aco,  the  radius  of  k(a  (o)  can  be 
shown  to  be: 

1 1/2 


radius  =  — - 

k(a  co)  2 


k(a  CO)  dco 


=  1 


1|k(co')||2 


(co’  -  l^(co’)  p  dco 


.1  _Ak 

1 


(4.15b) 


where 
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k(a(0)  1  2  =  I  k((o')  |  I2 


(4.15c) 


is  the  relationship  between  the  norms. 

From  the  above  analysis,  it  is  interesting  to  point  out  that  the  integral  reproducing  kernel 
window  transform  in  the  time  or  space  domain  [Eq.(3.5)]  and  in  the  frequency  or  wave  number 
domain  [Eq.(4.13)],  study  the  response  u  in  a  rectangular  time -frequency  or  space-wave  number 
window  given  by: 

[x  +  ax-aAk;x  +  ax  +  aAk]x(^-^;^  +  ^]  (4.16) 

The  above  window  narrows  to  pick  up  the  high-frequency  or  high  wave-number  phenomena  of  u 
and  widens  to  study  the  low_-frequency  or  low  wave-number  response.  This  suggests  that  we 
employ  a  flexible  time-frequency  or  space  wave-number  reproducing  kernel  window  to  define 
adaptive  refinement  of  the  local  response  of  u  around  any  point  x. 


4.4  Stability  analysis 


If  4>(y)  is  symmetric,  co  is  zero.  Therefore,  the  frequency  window  is  always  located  at 
(0  =  0,  and  its  frequency  band  becomes: 


(4.17) 


The  smaller  the  a,  the  larger  the  frequency  band;  the  larger  the  a,  the  smaller  the  frequency  band. 
This  implies  that  the  number  of  sampling  point  within  B(x)  must  satisfy  the  relation: 
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Ax  <  CON  X  7C  /  (^)  =  CON  X  2na  /  Ak 


(4.18) 


to  avoid  aliasing  (Brigham  [1974]).  Ax  is  called  the  sampling  rate.  The  constant,  CON,  will 
depend  on  the  so-called  frame  bounds  (Daubechies  [1992]): 

A<^|k(aco)p<B  (4.19) 

Ax 

The  positive  frame  bounds  coefficients  A  and  B  can  be  estimated  numerically.  If  A/B  is  close  to 
one,  CON  will  be  close  to  one.  Equation  (4.18)  is  referred  to  as  the  stability  condition  of  the 
reproducing  kernel  window  function  methods.  In  practice,  CON  is  chosen  much  less  than  one 
(Liu,  Zhang  and  Ramirez  [1991]). 

It  is  well  known  that  even  if  the  sampling  rate  satisfies  Eqs.  (4.18)  and  (4.19),  there  is  not 
a  good  choice  for  a  high  frequency  band  window  function  that  can  provide  accurate  frequency  and 
time  resolutions  of  u  simultaneously.  Therefore,  for  a  small  a,  an  intelligent  selection  of  the 
frequency  band  is  necessary  to  be  effective.  One  possible  way  to  employ  a  larger  sampling  rate  in 
high  frequency  or  high  wave  number  analysis  is  the  multiple-scale  method  proposed  by  Liu,  Zhang 
and  Ramirez  [1991].  In  this  approach,  the  response  is  divided  into  multiple  frequency  bands  via  a 
shifting  theorem.  If  ^>ab(’^)  ^  wavelet,  then  the  center  frequency  co  >  0.  By  breaking  a  up  into 

different  scales,  we  shall  have  a  multi-resolution  analysis  (Chui  [1992]). 

In  this  paper,  we  restrict  ourselves  to  the  analysis  of  a  single  frequency  band,  where 
^ab(*)  ^  chosen  to  be  the  scaling  function  so  that  ©  =  0.  Consequently,  depending  on  the  choice 
of  a  and  the  scaling  function,  d>(x),  the  integral  reproducing  kernel  window  transform,  (u,  kax),  is 
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a  measure  of  the  amount  of  change  of  u  at  the  location  x  +  ax  with  the  zoom-in  (smaller  a)  and 
zoom-out  (larger  a)  capability. 

4.5  Limitation  of  very  low  frequency/wave  number  analysis 

The  disadvantage  of  a  single  band  reproducing  kernel  analysis  is  that  the  method  will  break 
down  for  very  low  frequency/wave  number  analysis,  such  as  for  large  dilation  parameter  a.  This 
can  be  seen  by  interpreting  the  reproducing  kernel  Eq.(3.5)  as  a  continuous-time  or  space 
convolution: 


u‘>(x)  = 


{C(a.  X,  y)  a-‘  <I>(^^)}  u(y)  dy 


(4.20a) 


and  the  reproducing  kernel  is  identified  as: 

k(a,  y)  =  C(a,  y)  a  *  0(^)  (4.20b) 

U 


Applying  the  Fourier  transform  to  both  sides  of  Eq.(4.20a)  and  employing  the  Fubini  theorem 
(Daubechies  [1992])  to  the  right  hand  side  yields: 

u*'(o))  =  k(a,  ©)  u(co)  (4.21) 

In  order  for  u**(o))  to  approach  u(co),  k(a,  (o)  must  be  constructed  such  that  k(a,  (o)  =  1  for  all  O) ; 

ys. 

however,  this  will  violate  the  Rierr  'nn-Lebesgue  Lemma  that  k(a,  co)  =  0  when  to  approaches  ±oo. 
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An  approximation  of  the  convolution  identity  constructs  k(a,  to)  so  that  k(a,  to)  =  1  as  a  approaches 
to  zero.  It  is  well  known  that  the  family  of  Gaussian  functions: 

a-i  <I)(2.)  =  -L=  e  ,  a  >  0  (4.22) 

^  d.vz 

will  satisfy  the  above  conditions  provided  C(a,  x)  is  constructed  so  that  it  is  equal  to  one  when 
is  in  the  interior  of  the  domain  and  is  fairly  constant  when  moves  close  to  the  boundary.  With 

this  construction: 

k(a,  x)  =  a'*  <I>(^)  =  -7=  e  (4.23a) 

^  aVjt 


and  the  Fourier  transform  of  k(a,  x)  is: 


k(a,o))  =  e-“^“'/'^  (4.23b) 

With  the  Gaussian  function,  it  is  noted  that  a  is  identified  as  the  standard  deviation,  and  k(a,  CO) 
approaches  one  as  a  approaches  zero.  However,  when  a  is  large  Eq.(4.19)  will  break  down  since 
k(a,  co)  ^  1  unless  the  frequency  content  of  u(x)  is  close  to  zero,  which  is  very  restrictive. 

From  the  above  argument,  the  dilation  parameter  should  be  chosen  within  a  banded  range, 
say  amin  ^  a  ^  amax-  The  maximum  Umax  will  be  chosen  so  that  Eq.(4.23b)  is  close  to  one.  The 
minimum  amin  can  be  a  small  number  provided  the  stability  condition  Eq.(4.18)  is  met.  From  our 
numerical  experiments,  we  found  that  a  can  take  a  large  value.  We  believe  that  the  correction 
function  C(a,  x,  y)  indeed  improves  the  stability  as  well  as  the  accuracy  of  the  interpolation 
function.  Furthermore,  due  to  the  presence  of  C(a,  x,  y)  which  contains  the  (x  -  y)  term,  the 
reproducing  kernel  particle  interpolation  function  as  given  in  Eq.(3.7b)  will  further  increase  the 
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order  of  the  shape  function  by  one  so  that  the  convergence  rate  will  also  be  increased  by  one  when 
a  is  large!  This  wiU  be  discussed  further  in  the  next  section. 

We  shall  employ  the  Gaussian  function  and  the  cubic  spline,  which  is  a  good 
approximation  to  the  Gaussian  function,  as  the  window  function  <l>ax(y)  the  subsequent 

development  It  is  further  noted  that  both  Gaussian  and  cubic  spline  functions  are  scaling  functions 
used  to  generate  the  “Mexican  hat”  wavelet 


5.  Examples  of  Reproducing  Kernel  Window  Functions 


We  let  the  independent  functions  be: 


P(x)  =  {l,x} 

ID 

(5.1a) 

P(x)  =  {1,  Xi,  X2} 

2D 

(5.1b) 

P(x)  =  {1,  Xi,  X2,  X3} 

3D 

(5.1c) 

By  substituting  the  above  vector  of  independent  polynomials  into  Eqs.(3.3)  and  (3.6b),  the 
correction  function  C(a,  x,  y)  can  be  separated  into  two  terms; 

C(a,  X,  y)  =  Ci(a,  x)  +  C2(a,  (5.2) 

where  Cl  and  C2  =  {C21, ... ,  C2nsd}  ^  scalar  and  vector,  respectively,  which  can  be  defined 
in  terms  of  the  zero,  first  and  cross  moments.  It  is  noted  that  all  moments  are  in  general  a  function 
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of  the  location  of  the  window  function,  x,  eventhough  they  are  constants  if  B(x)  is  not  close  to  the 
boundary.  The  expressions  for  Ci  and  C2  in  ID,  2D,  and  3D  are  as  follows: 


One  dimension: 


Ci(a,  x)  = 


mn 

momii-m] 


(5.3a) 


C2(a,  x)  = 


mi 

momii-m] 


(5.3b) 


Two  dimensions: 


Ci(a,  x)  = 


_ mnm22-m^2 _ 

mo(m  Hm22-m]2)’(  m^m22-2m  im2m  i2+m|m  1 1) 


(5.4a) 


Cj,(a,x)=— -  ^  (5  4^, 

mo(mi im22-m-|2)-(  m-fm22-2m im2m i2+m|m 1 1) 


C22(a,  x)  = 


_ m2mii-mimi2 _ 

mo(miim22-m^2)'(  m]m22-2mim2mi2+m2mii) 


(5.4c) 


Three  dimensions: 


Ci(a,  x)  = 


2mi2m23m3i+mi  im22m33-miim2:^-m22m^i-m33m^2 

D 


(5.5a) 


C2i(a,  x)  = 


mi(m22m33-m^3)+m2(m23m3i-mi2m33)+m3(mi2m23-m3im22) 

D 


(5.5b) 
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C22(a,  x)  = 


m2(m33mi  i-m^  i)+m3(m3  im  i2-m23m  i  i)+m  i(m23m3i-m  121x133) 

D 


(5.5c) 


C23(a,  x)  = 


m3(mi  10122-01^2)+^  i(^i2ni23-ni3im22)+m2(ni3inii2-m23n^ii) 

D 


(5.5d) 


where  D,  the  Jacobian  of  the  4  x  4  M(x)  matrix,  is: 


D  =  mo(2mi2m23ni3i+mnm22m33-niiim^3-m22m|i-m33m^2) 
+mi(m^3-m22m33)+ni2(m3rm33mii)+m3(m^2*"ai  10122) 

+2mim2(m33mi2-m230i3i)+2m2m3(miim23-m3imi2) 

+2m3mi(m220i3rmi2m23) 

Because  of  the  special  properties  of  the  moments,  Ci(a,x)  =  1  aod  C2(a,x)  =  0  when  the 
window  function  is  within  R^;  whereas  Ci(a,x)  ^  1,  and  C2(a,x)  ^  0  when  the  window  function 

is  close  to  the  boundary.  It  is  particularly  important  to  point  out  that  if  a  is  large,  the  support  B(x) 
is  large.  Then  the  linear  term  appearing  in  Eq.(5.2),  C2(a,  x).(^  ^  ^),  plays  an  important  factor  in 

the  accuracy  and  convergence  rate  of  the  method. 

From  Eq.(3.5)  and  Eq.(3.6),  the  reproducing  kernel  function  and  the  approximation  uh(x) 
are  given  by: 


k(a,  X,  y)  =  C(a,  x,  y)  E(a) 


(5.6a) 


uh(x) 


k(a,  X,  y)  u(y)dRy 


(5.6b) 
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Therefore,  if  <h(r)  (r  =  1 1  x  -  y  1 12  )  is  the  cubic  spline  function,  the  order  of  polynominals  of  4>(r)  is 
equal  to  3.  In  a  Galerkin  formulation,  the  convergence  rates  of  the  solution  (L2  norm)  and  its  first 
derivatives  (HI  norm)  (see  numerical  examples  session  for  definitions)  are  expected  to  be  4  and  3, 

respectively.  However,  when  a  is  larger,  observing  from  the  discretized  Eq.(5.6a)  and  Eq.(5.2) 
(with  y  replaced  by  xj),  the  order  of  the  polynominal  of  the  reproducing  kernel  k(a,  x,  y)  can  be 

increased  by  one  so  that  the  L2  and  HI  convergence  rates  can  be  as  high  as  5  and  4,  respectively. 
This  unusual  phenomena  are  observed  in  our  ID  numerical  experiments.  The  convergence  rates  are 
much  higher  when  <I)(r)  is  replaced  by  a  Gaussian  function. 

Another  interesting  observation  can  also  be  abstracted  from  Eq.(5.2).  In  order  to  increase 
the  convergence  rate  by  one  order,  it  is  suggested  to  underintegrate  the  first  moments  so  that  it  is 
close  to  zero,  and  the  term  C2(a,  x).(^  ^  ^)  can  act  as  a  stabilization  term  to  the  SPH  methods.  At 

the  same  time,  this  stabilization  term  can  also  improve  the  accuracy  as  well  as  the  convergence  rate. 
One  way  to  achieve  this  is  to  use  Trapezoidal  Rule  to  integrate  the  M(x)  matrix  (that  is,  the 
moments)  at  each  x. 

For  higher  order  polynominals,  as  well  as  other  independent  functions,  P(x)  can  also  be 
similarly  investigated  for  this  class  of  RKPM  interpolation  functions;  however,  from  our  numerical 
experiments,  we  found  that  using  linear  polynominals  are  accurate  enough  for  most  purposes.  We 
align  found  that  linear  polynomials  give  numerically  more  stable  calculations  than  higher  order 
polynomials. 


6.  Similarities  among  SPH,  DEM,  EFGM  and  RKPM  Interpolation  Functions 

For  an  illustration  of  the  comparison  among  the  various  interpolation  funtions,  only  one 
dimensional  linear  polynomials,  P(x)  =  { 1,  x},  are  implemented.  The  window  function,  <I>(x),  is 
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chosen  to  be  positive  definite  and  an  even  function.  If  we  use  the  Trapezoidal  Rule  to  discretize  the 
reproducing  kernel  Eq.(3.5),  the  shape  function  N^(x)  of  the  RKPM  is  (a>0): 


N?(x)  =  j  Ci(x)  +  C2(x)  )  a-i  *1^^)  AMj  (6.1a) 


Ci(x) 


mil 


(6.1b) 


C2(x)^ - SiL—  (6.1c) 

momii-mf 

where  AMj  is  the  particle  mass.  If  mg  =  1 ,  m  u  0,  and  m  j  =  0,  the  SPH  interpolation  shape 
function  can  be  obtained.  As  a  matter  of  fact,  the  shape  functions  of  RKPM  and  SPH  are 
equivalent  in  an  interior,  but  there  is  a  big  difference  when  they  are  close  to  the  boundary.  Hence, 
the  SPH  methods  are  not  accurate  when  boundaries  are  present. 


Direct  differentiation  of  Eq.(6.1a)  gives: 


Nf..W  =  I  Cl'(x)  +  QXx)  ^ AM, 

+  (  C,(x)  +  C2(x)  AM, 


(6.2a) 


momn  -  m^ 


mu  (  mpmn  +  mom^^  -  2mimj ) 
(momn  -  m]  )2 


(6.2b) 


m 


momn  -  mt 


mi  (  rngmii  +  nrionijj  -  2mim.  ) 

(  momn  -  )2 


(6.2c) 
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In  the  above  expression,  a  superposed  " ' "  denotes  differentiation  with  respect  to  the  argument 

Without  going  into  details,  the  MLSM,  DEM  and  EFGM  interpolation  functions  can  be 
written  as  (Belytschko  et  al.[1993]): 

2 

Nj(x)  =  £  Pj(x)  [  A-i(x)  B(x)  ]j  j  (6.3a) 

j=i 

when  constant  and  linear  polynomials  Pj(x)  are  employed.  The  P(x),  A(x)  and  B(x)  matrices  are 
defined  by: 

pT(x)  =  {Pi,P2}  ;  Pi(x)  =  l  ;  P2(x)  =  x  (6.3b) 

A(x)  =  X  a-i<D(^)  P(xi)  pT(xi)  (6.3c) 

1=1 

B(x)  =  ja-itb(^)  P(xO,  ■  •  •,  a-^^)  P(xn))  (6.3d) 

where  n  is  the  number  of  points  in  the  neighborhood  of  x  for  which  the  weight  function 
^  0;  and  xj  are  the  nodal  coordinates  of  uj.  The  derivatives  of  the  EFGM,  Nj(x),  given 

by  Belytschko  et  al.[1993]  are  shown  as: 

Nf.xW  =  i  (Pj.x  +  Pj  [a;>B  +  A'1B,4,)  (6.4a) 

j  =  l 

and  the  derivative  of  the  A"^  matiix  is  given  by: 
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A;i  =  -A-‘A,xA-‘ 


(6.4b) 


The  derivative  of  Nj(x)  of  DEM  developed  by  Nayroles  et  al.[1992]  assumes  A  and  B 
constant  so  that: 


^.x(x)  =  t  Pj.x[A-‘B]jj  (6.5) 

j=i 

It  is  noted  that  no  particle  mass  AMj  nor  nodal  length  Axj  are  included  in  Eqs.(6.3)  through  (6.5). 
Furthermore,  the  A  and  B  matrices  need  to  be  computed  at  each  quadrature  point  x. 

Although  it  is  not  very  apparent,  an  interesting  result  arises  from  these  three  conditions. 

1)  The  nodal  coordinates, xp  are  equally  spaced. 

2)  The  Trapezoidal  Rule  is  used  to  numerically  integrate  M(x). 

3)  The  integration  weights  Axj,  AMj,  and  pj  in  Eqs.(6.1a)  and  (6.2)  are  all  set  equal  to  1. 

^  and  N?;*,  Eq.(6.1a)  and  Eq.(6.2),  can  be  shown  to  be  equivalent  to  Nj  and  Nj,x, 
Eq.(6.3a)  and  Eq.(6.4a),  respectively. 

With  the  above  assumptions,  Nj^x  can  be  similarly  defined  from  the  RKPM  interpolation 
functions  (with  AMj  =1): 
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(6.6) 


N?.,(x) - ,  ^  'iMj 

almomii-mfj 


am, 

aimomii-m^J 


Comparing  Eq.(6.2)  with  Eq.(6.6),  depending  on  d)(x),  the  derivative  N?,*  might  not  be  an 
accurate  approximation  to  N^x  >  especially  for  large  a  values.  For  example,  let  us  use  the  cubic 
spline  functions: 


+  for  0<-^<^ 


2  Ax2  Ax^ 


(6.7a) 


.  Ax 


<I,(r)  =  i.  4l  + ill.  fori<-E-<l 

^  Ax  Ax2  3Ax^  ^  Ax 


(6.7b) 


r  =  ll  y  -  xlb 


(6.7c) 


as  the  window  function,  and  the  dilation  parameter  is  set  such  that  a  =  for  j>0.  The 

parameter  j  =  0  is  adjusted  so  that  Ax  is  right  at  the  stability  limit  (see  Eq.(4.18)  and  Eq.(4.19)). 
This  corresponds  to  the  smallest  time  bandwidth  or  the  largest  frequency  bandwidth  of  the  window 
function.  Let  us  consider  a  set  of  21  equally  spaeed  nodes  with  Ax  =  0.3  representing  the  domain, 
0  ^  X  ^  6.0.  If  we  use  a  trapezoidal  rule  to  discretize  the  domain,  then  Axj  =  Ax2i  =  0.15,  and 
Ax2  =  Ax3  =...=  Ax2o  =  0.3.  With  this  discretization  and  p  =  1,  the  shape  function  given  in 
Eq.(6.1a),  its  derivatives  using  Eq.(6.2)  [exact]  and  Eq.(6.6)  [approximation]  are  depicted  in 
Figures.!  and  2,  respectively  for  nodes  1,  10,  and  21.  As  might  be  seen  for  j  =  0  (right  at  the 
stability  limit),  the  shape  functions  and  its  derivatives  look  similar  to  those  for  the  usual  linear  finite 
elements;  however,  the  derivatives  of  RKPM  shape  functions  are  continuous  and  they  try  to 
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reproduce  the  finite  element  discontinuous  derivatives.  For  j  =  1.0  and  j  =  2.0,  there  is  not  much 
difference  between  the  two  derivatives  especially  for  node  10  where  the  support  is  within  the 
domain.  There  is  a  large  difference  between  the  two  derivatives  when  the  shape  functions  are 
located  at  the  boundaries  (x  =  0.0  and  x  =  6.0).  This  difference  will  produce  inaccurate  derivatives 
in  Eq.(6.6),  and  the  solution  deteriorates. 


If  we  pick  <X>(r)  as  the  Gaussian  function  such  that: 


-1  ^  y  .)  =  ^_e-(x-y)-/a-  ,  a>0 


(6.8a) 


a=  j  ^0  and  = 

afiz 


(6.8b) 


With  O  as  the  standard  deviation,  the  exact  derivative  given  in  Eq.(6.2)  and  the  approximate 
derivative  given  in  Eq.(6.6)  become,  respectively: 


Nf,,(x)  =  j  Ci'(x)  +  CiXx)  ^  AMj 

-  \  (  C.(x)  ^  C,(x)  AM, 


(6.9a) 


ns;,(x)  =  -  ,  AM. 

avmom[i-mjJ 


mi 


momi 


i-m]) 


a-i<D(^)  AMj 


(6.9b) 
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When  the  derivatives  of  the  shape  function  are  not  evaluated  close  to  the  boundaries  (mi=0),  the 
two  derivatives  are  represented  with  the  same  function  but  different  coefficients. 

When  the  derivatives  are  evaluated  close  to  the  boundaries,  the  two  derivatives  are  again 
represented  by  the  same  Gaussian  function,  and  from  numerical  experiments,  the  quadratic  term, 
»  appearing  in  Eq.(6.9a)  does  not  play  an  important  role.  Therefore,  the  approximation  in 

Eq.(6.9b)  is  a  good  approximation  to  Eq.(6.9a).  This  is  further  elaborated  in  the  next  section. 


7.  Numerical  Experiments 


We  employ  a  one  dimensional  Laplacian  type  equation; 


U.XX  +  2  s2  sech^[s  (x  -  3)]  tanh[s  (x  -  3)]  =  0 

(7.1a) 

with  essential  boundary  conditions: 

u(0)  =  -  tanh(3s) 

(7.1b) 

u(6)  =  tanh(3s) 

(7.1c) 

The  parameter  s  controls  the  degree  of  localization  of  the  gradient  of  u  (u,, 

has  an  increasing  gradient.  The  exact  solution  is; 

.).  As  s  increases 

u(x)  =  tanh[s  (x  -  3)] 

(7. Id) 
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In  a  numerical  approximation,  we  employ  a  Galerkin  formulation  of  Eq.(7.1a)  and  the  boundary 
conditions  are  enforced  via  the  standard  Lagrange  multiplier  approach.  For  a  detail  description  of 
this  problem  and  the  implementation  of  the  boundoi'y  conditions  (Belytschko  et  al.  [1993]). 

We  shall  utilize  the  cubic  spline  and  Gaussian  function  described  in  section  6  as  window 
functions.  For  simplicity,  linear  polynomials,  P(x)  =  { 1,  x},  and  s  =  10  are  used  throughout.  Five 
different  discretizations  with  Ax  =  0.3  (21  nodes),  0.15  (41  nodes),  0.075  (81  nodes),  0.0375 
(161  nodes),  and  0.025  (241  nodes)  are  solved.  It  is  noted  that  similar  to  SPH,  DEM,  and  EFGM, 
RKPM  does  not  require  an  element  nor  element  connectivities.  The  standard  L2  and  HI  error 
norms  are  defmed  as: 


(L2  norm)^  = 


(HI  norm)^  = 


(7.2a) 


(7.2b) 


The  rates  of  convergence  are  defined  as  the  slopes  p  and  q  appearing  in  the  In  (error)  vs.  In  (Ax) 
equations: 


ln(L2  noiTn)  =  InGi  +  p  InAx 

(7.3a) 

ln(Hl  norm)  =  lnG2  +  q  InAx 

(7.3b) 

The  smaller  the  constants  Gj  and  G2,  the  more  accurate  the  method.  Also,  the  higher  p  and  q,  the 
faster  the  rate  of  convergence.  To  obtain  the  convergence  plots,  ten  and  twelve  Guass  quadrature 
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points  are  used  to  integrate  the  matrices  and  the  errors,  respectively.  Nevertheless,  only  three  to 
five  Guass  quadrature  points  are  sufficient  to  integrate  the  matrices  accurately.  The  Trapezoidal 
Rule  is  used  to  integrate  all  the  moments  and  the  M(x)  matrices. 

The  exact  solution  of  u,x  and  the  SPH  Galerkin  approximation  using  a  Gaussian  window, 
81  nodes  (Ax  =  0.075),  o  =  Ax/CO.  IVtT)  ,  and  Ax/fO.SVrT)  is  depicted  in  Figure  3.  It  shows  the 
SPH  solution  depends  very  much  on  the  dilation  parameter  (or  the  standard  deviation).  As  a  matter 
of  fact,  for  a  given  Ax,  the  larger  the  o  (the  larger  window),  the  better  the  gradient  is  at  the  center. 
However,  both  choices  of  o  give  bad  approximations  of  the  gradient  of  u  at  the  boundaries.  This 
confirms  that  SPH  interpolation  functions  do  not  work  well  with  boundaries. 

The  L2  and  HI  norms  plots  for  the  spline  window  function  are  in  Figures  4  and  5, 
respectively.  As  can  be  seen  in  Figures  4b  and  4d,  and  Figures  5b  and  5d,  Eq.(6.2)  works  better 
than  Eq.(6.6)  when  a  is  large  (i.e.,  j  is  large).  Similar  conclusions  can  also  be  drawn  from  Figures 
4a  and  4c,  and  Figures  5a  and  5c.  One  interesting  observation  is  that  the  L2  and  HI  convergence 
rates  (see  Figures  4c  and  5c)  can  be  as  high  as  5.7  and  3.9,  respectively.  This  confirms  our 
analysis  in  Section  6. 

Similarly,  the  L2  and  HI  convergence  rates  plots  shown  in  Figures  6  and  7  are  produced 
using  the  Gaussian  window  function.  The  standard  deviation  and  the  dilation  parameter  are  chosen 
such  that: 


and  a  =  2^^^’^^cy  (7.4) 

oViT  VI 

As  predicted  in  Section  6,  the  two  formulas  for  the  derivatives  give  virtually  identical  results.  It  is 
interesting  to  point  out  that  the  L2  and  HI  convergence  rates  can  go  as  high  as  15.92,  and  14.91, 
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respectively.  Finally,  in  Table  1,  the  peak  values  of  are  compared  with  the  standard  linear  finite 
element  method  (FEM)  and  the  exact  solution.  The  RKPM  is  able  to  capture  the  high  resolution  of 
the  steep  localized  gradient  using  a  flexible  space-wave-number  Gaussian  window  function. 


8.  Conclusions 

In  this  paper,  a  review  of  the  mesh  or  grid  free  interpolation  functions  is  presented.  The 
dilation  and  translation  of  a  window  function,  the  integral  window  transform,  and  the  SPH 
interpolation  kernel  functions  are  reviewed.  By  understanding  the  merits  and  deficiency  of  the 
SPH,  MLSM,  DEM,  and  EFGM  methods,  a  new  continuous  reproducing  kernel  particle 
interpolation  kernel  function  is  derived  in  terms  of  a  fiexible  time-frequency  or  a  space-wave 
number  localized  window  function.  Comparing  to  the  SPH  interpolation  function,  a  continuous 
correction  function  to  the  SPH  methods,  which  is  composed  of  the  various  moments  of  the 
window  function,  is  identified. 

The  effect  of  the  dilation  parameter  and  the  stability  condition  of  this  new  discretized 
reproducing  kernel  particle  interpolation  function  are  discussed.  The  convergence  rate  of  this  class 
of  RKPM  interpolation  functions  using  a  Galerkin  method  is  shown  to  be  at  least  one  order  higher 
that  that  of  the  window  function  when  the  dilation  parameter  is  large.  Since  the  correction  function 
and  the  window  function  can  be  chosen  to  be  smooth,  the  solution  as  well  as  its  derivatives  are 
continuous  throughout  the  entire  domain  of  interest,  unlike  the  usual  finite  element  methods.  The 
numerical  experiments  confirm  the  theoretical  analysis  presented  in  this  paper. 
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Table  1.  Comparison  of  derivative  peak  values  among  RKPM,  FEM,  and  exact  solution. 


Fig.  1.  Plots 


Derivatives  of  shape  functions 


Fig.  2.  Derivatives  of  interpolation  functions  using  exact  and  approximate  formulas  at 
node  1  (x  =  0.0),  node  10  (x  =  2.7),  and  node  21  (x  =  6.0). 


81  nodes 


Spline  (Approximation )  Spline  ( Exact ) 


Fig.  4.  L2  norm  and  convergence  plots  using  approximate  and  exact  formulas  for  the  derivatives 
of  the  interpolation  function  using  a  (cubic  spline)  window  function. 


Spline  (Approximation )  Spline  ( Exact ) 


Hi  norm  and  convergence  plots  using  approximate  and  exact  formulas  for  the  derivatives 
of  the  Interpolation  function  using  a  (cubic  spline)  window  function. 
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Fig.  6.  L2  norm  and  convergence  plots  using  approximate  and  exact  formulas  for  the  derivatives 
of  the  interpolation  function  using  a  (Gaussian)  window  function. 
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Abstract 


This  analysis  explores  a  Reproducing  Kernel  Particle  Methods  which  incorporates  several 
inviting  features.  The  emphasis  is  away  from  classical  mesh  generated  elemen^  in  favor  of  a 
mesh  free  system  which  only  requires  a  set  of  nodes  or  particles  in  space.  Using  a  Gaussian 
distribution,  flexible  window  functions  are  implemented  to  provide  refinement  in  the  solution 
process.  It  also  creates  the  ability  to  analyze  a  specific  frequency  range  in  dynamic  problems 
reducing  the  computer  time  required.  This  advantage  is  achieved  through  an  increase  in  the 
critical  time  step  when  the  frequency  range  is  low  and  a  large  window  is  used.  The  stability  of 
the  window  function  is  investigated  to  provide  insight  on  Reproducing  Kernel  Particle  Methods. 
Furthermore,  there  are  no  explicit  elements  in  the  formulation,  allowing  the  derivatives  to  also  be 
continuous,  or  C°°.  The  analytic  theory  is  confirmed  through  numerical  experiments  by 
performing  reconstructions  and  solving  an  elastic-dynamic  one  dimensional  problem. 

1.  Introduction 

There  is  always  a  drive  to  find  new,  more  advantageous  ways  to  analyze  problems  using 
numerical  methods.  Typical  finite  elements  use  linear  or  quadratic  shape  functions  to  define  the 
response  within  each  element.  For  large  deformation  or  high  frequency  problems,  the  elements 
must  be  very  small  to  accurately  predict  the  characteristic  response.  To  avoid  these  problems,  the 
Reproducing  Kernel  Particle  Method  (RKPM)  advantages  are  exploited. 

There  is  no  explicit  mesh,  so  mesh  creation  time  is  saved.  Since  a  mesh  is  not  required, 
there  are  not  any  problems  due  to  mesh  entanglements  allowing  for  large  deformations  pd 
unrestrained  movement  of  nodes.  Not  only  is  mesh  creation  time  saved,  but  mesh  recreation  time 
is  eliminated  since  to  refine  the  problem  in  an  area  of  interest  one  needs  only  to  add  points  in  the 
interesting  region.  The  method  for  meshless  finite  elements  based  on  the  moving  least  square 
interpolant  was  used  by  Lancaster  and  Salkauskas  (1981)  for  graphics,  but  has  been  given  new 
life  recently  by  Nayroles  et  al.  (1992)  with  application  to  Diffuse  Element  Methods.  This  was 
followed  by  the  Element  Free  Galerkin  Method  (EFGM)  by  Belytschko  et  al.  (1993)  that 
improved  upon  Nayroles  method,  including  a  correction  to  the  derivative  of  the  shape  function. 
There  has  also  been  work  performed  on  Smooth  Particle  Hydrodynamic  (SPH)  method.  The  SPH 
method,  developed  by  Gingold  and  Monaghan  (1977)  and  others,  provides  a  meshfree 
environment  but  it  has  some  difficulty  creating  accurate  solutions  on  the  boundaries  or  when  a 
small  number  of  particles  is  used  [Monahgan  (1988)].  The  SPH  method  is  similar  in  basic 
construction  as  the  Gaussian  Reproducing  Kernel  (GRK)  to  be  presented  here,  but  it  lacks  key 
features  of  the  GRK. 

In  this  Gaussian  Reproducing  Kernel,  the  procedure  used  is  based  on  continuous 
reproducing  kernel  particle  construction.  It  can  also  be  viewed  as  a  continuous  Least  Square 
polynomial.  However,  by  exploiting  the  moment  definitions  of  the  flexible  window  function,  the 
reproducing  kernel  can  be  reduced  to  a  simplified  form  so  its  properties  can  be  investigated. 
Second,  a  localized  flexible  window  function  is  incorporated  by  translating  the  function  across 
the  entire  domain  to  reproduce  the  response.  The  window  function,  part  of  the  shape  function,  is 
controlled  by  two  parameters.  Unlike  the  typical  finite  elements  which  only  have  one,  the  two 
parameters  allow  a  greater  problem  solving  ability.  The  flexible  window  allows  for  response 
frequencies  or  wave  numbers  to  be  selectively  reproduced  in  the  numerical  approximation. 
Essential  to  the  development  of  this  method  is  an  understanding  of  the  stability  limits  of  the 
flexible  window  function,  as  well  as  the  critical  time  step  in  dynamic  analyses.  The 
aforementioned  properties  of  the  flexible  window  function  and  reproducing  kernel  create 
interesting  characteristics  for  the  kernel  stability.  Determining  equations  are  developed  for  both 
of  these  stability  limits. 

The  reproducing  kernel  in  this  derivation  turns  out  to  be  similar  to  a  free  Lagrangian 
particle  method  [Libersky  (1990)]  with  one  major  difference-the  development  of  a  correction 
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function.  The  primary  motivation  behind  the  correction  function  is  to  provide  accurate  solutions 
at  the  boundaries,  but  by  careful  integration  techniques  it  is  also  possible  for  the  correction 
function  to  provide  stability  to  the  solution.  Free  Lagrangian  methods  provide  accurate  solutions 
in  the  interior  of  the  problem  when  the  number  of  particles  is  large,  but  they  do  not  provide  a 
means  to  get  an  accurate  solution  near  the  boundary.  This  method  incorporates  correction 
functions  that  are  relatively  dormant  in  the  interior  and  then  provide  correction  on  demand  at  the 
boundaries. 

Through  the  implementation  of  a  window  function  and  the  knowledge  of  the  Fourier 
transform,  it  is  possible  to  develop  a  new  type  of  shape  function  that  can  still  be  used  in  the  usual 
Galerkin  formulation.  The  derivative  of  the  shape  function,  and  thus  reproducing  kernel,  can  be 
obtained  by  direct  differentiation.  The  development  of  the  proposed  shape  function  will  be 
derived  in  detail  later,  but  for  now  we  describe  its  characteristics.  The  two  parameters  in  the 
shape  function  provide  the  ability  to  translate  and  dilate  the  window  function.  Translation  is 
required  to  move  the  window  function  around  the  domain  since  the  window  functions  themselves 
have  a  compact  support.  The  ability  to  translate  replaces  the  need  to  define  elements.  The 
dilation  parameter  is  used  to  provide  refinement  to  the  solution  by  reducing  the  number  of 
calculations  necessary  to  find  a  solution.  The  larger  the  dilation  parameter,  the  smaller  the 
frequency  band  is  in  the  solution,  and  the  larger  the  critical  time  step  becomes  in  dynamic 
analyses.  The  refinement  parameter  transformation  between  the  time  and  frequency  domain  (or 
space  and  wave  number)  controls  the  solution  space.  This  introduces  the  ability  to  choose  the 
size  of  the  frequency  or  wave  number  range  in  the  calculation. 

2.0  Development  of  the  Reproducing  Kernel 

Prior  examples  of  reproducing  kernels  can  be  found.  The  most  obvious  may  be  its  use  in 
the  Fourier  Transform  Analysis.  The  Fourier  Transform  motivates  this  analysis  since  this  is 
where  the  concept  to  analyze  specific  frequency  bands  incorporated  into  this  type  of  reproducing 
kernel  analysis  originates. 

The  general  form  of  a  reproducing  kernel  is  a  class  of  functions  that  output  the  function 
itself  when  the  integral  over  the  domain  is  performed.  The  Fourier  Transform  is  an  excellent 
example  of  a  reproducing  kernel.  Its  intricacies  are  defined  by  the  following  set  of  equations 

where  the  Fourier  Transform,  f{co)  .  of  a  function,  f{x),  is  defined  by  equation  (2.0.1)  and  the 
inverse  Fourier  Transform  is  described  by  (2.0.2). 


fix) 


fix)  = 


e-ix(of(^x)dx 


gixmj^x)  da 


(2.0.1) 


(2.0.2) 


The  bounds  of  the  Transform  can  be  easily  described  if  the  interesting  spectrum  contains  a  single 
frequency  band.  The  bounds  of  the  Fourier  transform,  f{a),  are  given  in  equation  (2.0.3). 
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It  is  now  possible  to  show  that  by  performing  the  inverse  transformation  on  the  transformation  to 
obtain  the  original  function,  the  following  solution  is  obtained  where  fix)^0  only  in  the 
domain  x  e  [0,  L] . 
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This  representation  of  a  reproducing  kernel  using  a  sine  function  is  presented  in  Chui  (1992)  as 
well  as  Liu  (1991).  In  Liu’s  analysis,  the  multiple  scale  analysis  predicts  a  critical  time  step 
according  to  the  following  equation. 


At<^  (2.0.5) 

Aa 

This  is  used  as  a  guideline  since  equation  (2.0.5)  determines  the  sampling  rate  to  prevent  alining 
using  the  reproducing  kernel  for  a  given  frequency  band,  Dw.  It  is  believed  that  after  discretizing 
the  continuous  reproducing  kernel,  a  similar  form  exists  for  the  stability  condition  of  the 
reproducing  kernel  in  space  (see  section  4). 

2.1  Window  Function  Selection 

The  ideal  window  function,  F(x),  is  chosen  such  that  the  following  two  conditions  are 
satisfied.  Its  integral  over  the  domain,  Rx,  must  be  unity,  and  it  must  also  be  orthogonal. 


1)  <P{x)dR,=\ 

JRx 


(2.1.1) 


over  the  entire  domain,  rather  the  value  is  zero.  Furthermore,  most  wavelets  cannot  capture 
constant  or  linear  terms,  so  they  cannot  reproduce  the  simplest  of  functions  very  well.  For  these 
reasons  wavelets  are  not  considered  here,  but  appealing  characteristics  are  used  from  the  wavelet 
analysis,  Chui  (1992).  Splines  are  always  an  option  since  we  can  custom  tailor  them,  but  they 
produce  unnecessary  irregularities.  The  Gaussian  function  is  not  orthogonal,  but  the  Gaussian 
function  is  used  in  this  analysis  because  of  its  special  properties.  It  is  also  important  to  keep  in 
mind  that  the  integration  over  the  entire  domain  will  be  less  than  one  on  the  boundaries. 

2.2  Development  of  Window  Function 

Like  the  Fourier  series,  we  want  to  be  able  to  reconstruct  any  function  or  response  u(x)  by 
a  series  of  window  functions,  F(x).  Since  F(x)  is  a  localized  function,  it  becomes  necessary  to 
translate  the  function  to  represent  the  entire  response.  This  is  performed  by  inserting  the 
argument,  s  -x,  in  the  function.  Now  the  response  is  represented  by  the  following  integral. 


JC  (!>(s  -x)dR:,=  1 

Rx 


(2.2.1) 


Where  C  is  the  correction  function  to  be  defined  later  in  the  construction,  section  2.4.  Or  it  could 
be  an  arbitrary  constant  as  in  free  Lagrangian  particle  methods,  Libersky  (1990). 

To  the  window  function  argument,  it  is  also  necessary  to  add  the  dilation  or  refinement 
parameter,  r.  This  is  incorporated  by  dividing  the  window  function  argument  by  this  refinement 
parameter.  For  notational  purposes  this  definition  is  used  for  the  window  function: 


C'P„(x)dR.=  \  C0(x)dR,  =  \  (2.2.2) 

JRx  Jrx 

An  intuitive  sense  of  the  refinement  parameter  can  be  revealed  by  thinking  of  the  parameter  as 
the  standard  deviation,  when  the  Gaussian  distribution  is  used  for  the  window  function.  The 
additional  constant,  r  ■^,  is  the  proper  constant  only  for  one  dimension,  Liu  (1993).  This  constant 
scales  the  window  function  such  that  integral  over  the  domain  of  the  window  equals  one 
according  to  equation  (2.2.2).  It  also  is  useful  to  derive  the  moment  equations  in  section  2.4.1. 


C  l0(^)dRx  =  l  (2.2.3) 

Jrx 

This  completes  the  development  of  the  window  function  to  be  used  in  this  1-D  analysis. 

2.3  The  Reconstruction  Equation 

The  concept  that  any  function  can  be  represented  as  a  sum  of  linearly  independent 
functions  initiates  the  analysis,  starting  with  the  following  definition. 

u{x)  =  P{x)d  ^  (2.3.1) 

Where  P(x)  =  [P](x),  P2(x), ...  ,  Pn(x)]  which  is  any  number  of  linearly  independent  functions 


dimensional  case,  P(x)  =  [1,  x],  or  for  a  quadratic  one  dimensional  case,  P(x)  =  [1,  x,  x^],  etc. 
This  Moving  Least  Square  Interpolant  type  of  reconstruction  has  been  used  by  Nayroles  (1992) 
and  Belytschko  et  al.  (1993),  for  the  Element  Free  Galerkin  Method. 

It  is  possible  to  solve  for  the  unknown  coefficients  d  by  using  the  window  function.  The 

variable  x  is  changed  to  s  in  equation  (2.3.1),  and  then  both  sides  are  premultiplied  by  r  P^is) 
and  the  integral  window  transform  is  applied,  Liu  (1993])  (Note:  r  is  the  density  and  it  is 
included  here  to  make  the  transition  to  a  larged  deformation  case  easier.)  The  integral  window 
transform  multiplies  both  sides  by  the  window  function,  F rx>  then  integrates  over  the 
domain. 


JpP^(s)u(s)  <f>rsds 

Rx 
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(2.3.2) 


c(x)=|  pP'^{s)P{s)  d>rxds  (2.3.3) 

JRx 

Equation  (2.3.3)  is  used  in  the  next  section  to  complete  the, construction  of  the  correction 
function.  (In  that  section  the  merits  of  C  are  discussed  which  shows  why  the  term  is  sometimes 
referred  to  as  the  boundary  correction  function.)  By  making  this  definition  for  C,  the  solution  for 
dl  can  be  substituted  back  into  (2.3.1),  obtaining  the  reconstruction  equation  in  (2.3.5)  which  can 
also  be  written  as  one  integral  (2.3.6). 


=  rkx) 

JRr 


pP^(s)u(s)  d>rxds 


(2.3.4) 


'\x)  I 

Jr 


uix)  =  Pix)  C^{x)  I  p  P^is  )  u{s  )  0rx  ds 


(2.3.5) 


m(x)=|  P  P(x)  C^ix)  P^(s)  uis)^  <Pl^-^)ds 
Jrx 


(2.3.6) 


2.4  Correction  Function 

Comparing  the  reconstruction  equation  in  (2.3.6)  to  the  SPH  method  establishes  that  the 

only  difference  is  the  appearance  of  the  P(x)  C^{x)  P^s)  term  in  the  GRK  method.  This  term  is 
defined  as  the  correction  function,  and  its  merits  are  analyzed  in  this  section.  For  simplicity,  the 
characteristics  of  this  function  are  derived  here  in  one  dimension  for  linear  polynomials  P(x),  but 
it  can  be  shown  to  be  valid  for  multiple  dimensions  as  well,  Liu  (1993).  Expanding  equation 
(2.3.3)  reveals  the  following.  r 


The  inverse  of  this  matrix  will  be  computed  along  with  the  r,  P(x)  and  P'^(s)  from  the 
reconstruction  formula  in  equation  (2.3.6)  to  form  the  correction  function.  The  entire  term, 

P(jc)  C^(x)  P^is),  function  simplifies  to  one  number  regardless  of  the  number  of  terms  used  for 
P(x)  and  P'^(s) . 


i 


mo  =  \  p^iz)  dz 
Ibm 


mi=  \  pzd>{z)dz 


m\i  =  1  pz'^  d>(z)  dz 


(2.4.1. 1) 


(2.4. 1.2) 


(2.4. 1.3) 


These  moment  equations  are  integrated  over  the  region  B(x),  where  B(x)  is  the  region 
where  the  window  function  is  non-zero.  The  calculation  could  be  performed  over  the  entire 
domain;  however,  there  are  many  unnecessary  calculations  involved.  In  order  to  determine  the 
support  B(x)  using  the  Gaussian  function,  a  three  standard  deviation  criterion  can  be  enforced 
which  insures  the  calculation  of  99.7%  of  the  total  area.  The  moment  equations  have  the 
following  characteristics  that  provide  for  an  accurate  solution  near  the  boundary,  and  lie 
relatively  dormant  in  the  interior. 


Interior  Region 

Near  Boundaries 

mo 

=  1 

<  1 

mi 

=  0 

mil 

=  r^ 

^  r^ 

Although  our  initial  indication  from  the  correction  function  was  that  it  vanished  in  the 
interior,  it  is  now  known  that  through  careful  integration  of  the  function  and  its  inverse  the 
function  can  have  a  profound  effect  on  the  stability  of  the  kernel,  see  figures  4  and  5.  It  is  noted 
that  the  stabilization  effect  is  much  more  pronounced  in  the  data  when  the  number  of  nodes  is 
relatively  small. 

2.4.2  Final  Correction  Function  Form 

The  solution  to  P{x)  C^ix)  P^s)  can  be  written  in  the  simplified  form  shown  below 
through  manipulation.  Equation  (2.4.1)  can  be  inverted  and  substituted  back  into  the 
reconstruction  in  equation  (2.3.6)  to  reveal  a  continuous  reproducing  kernel  for  the  function  u. 
This  is  the  final  reconstruction  equation  to  be  used  in  this  analysis!  (Note  that  the  correction  term 
is  simplified  into  the  sum  of  two  terms  which  are  defined  in  the  following  manner.) 


Ci  = 


^11 

(momii  -  mf) 


(2.4.2.1) 


mi 


^2=7 
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(2.4.2.2) 


(2.4.2.3) 


«(l)  =  [Cl  +  C2(i^)  ]  u(s)  j:  {^)  ds 
JRx 

This  is  the  reconstruction  equation  which  is  the  basic  building  block  of  this  method,  and 
its  intricacies  are  revealed  throughout  this  investigation.  In  order  to  proceed,  it  is  n^essary  to 
define  the  discretized  form  in  order  to  implement  the  reconstruction  into  the  numeric  solution 
process  later. 


NP 

uHx)  =  X  [Ci  +  C2(^^)  ]}  i4xj)  AMj  (2.4.2.4) 

J=\ 

The  discretized  reconstruction  equation  can  be  written  in  a  more  familiar  form  to  perform 
the  numerical  analysis.  The  shape  function  is  defined  as  Nj  in  equation  (2.4.2. 6). 


NP 

uKx)  =  X 


y=i 


(2.4.2.5) 


Nj{x)  =[  Cxix)  +  Ciix)  (^^)]  r-i  AMj 


(2.4.2.6) 


By  analyzing  the  properties  of  the  moments,  in  a  continuous  case  the  values  for  Ci  and 
Cl  are  found  to  equal  1  and  0  respectively  in  the  interior  of  a  solution  (assuming  a  sufficiently 
large  number  of  particles),  but  definitely  not  equal  to  these  values  on  the  boundary.  In  this 
continuous  integral  case,  the  GRK  will  be  identical  to  the  SPH  in  the  only  in  the  interior. 
However,  in  the  discretized  form,  the  effect  of  the  correction  function  will  depend  on  the 
integration  technique  used  as  well  as  the  number  of  particles  in  the  integration  domain.  These 
variables  will  determine  whether  the  correction  function  will  only  act  on  the  boundaries  or 
enhance  the  stability  of  the  solution  throughout  the  entire  domain. 


2.5  Gaussian  Reproducing  Kernel  Formulation 

It  is  important  at  this  point  to  define  the  Gaussian  function  used  in  this  analysis  at  this 
time.  The  refinement  parameter  is  defined  to  contain  a  measure  of  normalization  so  that  a  given 
dilation  of  the  window  function  always  contains  the  same  number  of  nodes,  regardless  of  particle 
density.  The  refinement  parameter  can  also  be  recognized  as  the  standard  deviation  in  the 
Gaussian  equation. 


J-  -1  e  -  Cc  - 

(2.5.1) 

(2.5.2) 

Defining  the  Gaussian  function  in  this  way  has  the  advantage  of  maintaining  the  same  number  of 
nodes  for  support  while  changing  the  distance  between  nodes,  but  the  window  and  shape 
functions  will  change  while  changing  nodal  coordinates  as  well  as  the  j  refinement.  The 
parameter  j  can  be  any  real  number  but  stability  limits  are  set  in  section  5.  In  prior  analyses,  the 
following  definition  was  found  to  be  optimal  ^Liu  (1993)],  so  it  is  chosen  as  a  starting  point  for 
this  analysis. 


2.6  Dynamic  Frequency  Analysis 

One  of  the  interesting  aspects  of  this  method  is  its  ability  to  preconceive  the  frequency 
range  studied  in  the  analysis.  The  frequency  range  captured  can  be  determined  by  calculating  the 
Fourier  Transform  of  the  window  function.  The  window  function’s  shape  can  be  changed  by  the 
refinement  parameter,  r,  allowing  for  adaptability  in  the  solution  process.  Even  though  the 
Fourier  Transform  of  the  Gaussian  function  is  also  bell  shaped,  a  cutoff  frequency  for  the  banded 
window  has  been  found  to  be  1/r ,  Liu  (1993).  This  method  matches  the  area  under  the  Gaussian 
function  with  the  area  of  the  straight  box  created  by  the  cutoff  frequency.  The  approximation,  1/r, 
for  the  highest  accurately  reconstructed  frequency  inside  the  window  function  enables  an 
understanding  of  the  limitations  of  this  method  in  its  straight  form. 

Unfortunately  for  this  undeveloped  implementation  of  this  method,  the  frequency 
window  is  always  centered  around  w  =  0.  This  means  that  we  cannot  selectively  consider  only 
higher  frequency  bands;  it  is  necessary  to  capture  all  the  frequencies  below  the  highest  frequency 
of  interest.  The  ability  of  this  method  to  capture  high  frequencies  is  also  limited  by  the  stability 
of  the  reproducing  kernel  itself.  Although  the  number  of  nodes  in  any  one  window  function  is 
variable,  there  must  be  at  least  two  nodes  to  remain  stable.  This  is  necessary  in  order  to  have 
connectivity  between  window  functions.  Without  the  variable  connectivity  being  always  greater 
that  two  nodes,  the  response  cannot  be  translated  along  to  the  adjacent  nodes. 

It  is  hoped  that  the  Multiple  Scale  Methods  by  Liu  (1991),  can  be  included  in  subsequent 
analyses  to  shift  the  interested  frequency  bands  away  from  the  origin,  removing  the  unnecessary 
calculation  of  the  frequencies  in  between  the  interesting  areas.  Another  approach  is  to  use 
wavelets  to  capture  the  high  frequency  bands  which  is  being  investigated  by  Liu,  Oberste- 
Brandenburg,  and  Chen. 

3.0  Given  Function  Reconstruction 

The  reproducing  kernel  is  used  to  perform  given  function  reconstruction,  or  curve  fitting, 
as  a  demonstration  and  evaluation  technique  before  it  is  implemented  into  a  mesh  free  Galerkin 
type  formulation.  This  is  done  by  using  the  discretized  form  of  the  reconstruction  equation 
(2.3.6)  to  reproduce  a  known  function.  The  Trapezoidal  Rule  is  used  here  for  integration  despite 
the  fact  that  it  is  known  to  under  integrate.  The  dilation  parameter  has  a  profound  effect  on  the 
shape  function.  Examples  of  shape  functions  for  several  values  of  the  dilation  parameter  are 
shown  below.  Note  that  for  j  =  -2  the  window  function  approaches  the  ordinary  finite  element 
shape  function. 
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Figure  1  Shape  function  dilations  for  j  =  -2, 0,  1,2 


Figure  2  Reconstruction  of  sin(a)x)  where  co  =  1  with  21  nodes,  demonstrating  the  advantage 

of  the  GRK  over  SPH  near  the  boundaries. 


Figure  1  portrays  that  as  the  dilation  parameter  decreases  in  magnitude,  the  window 
function  approaches  a  Dirac  delta  function.  If  it  was  possible  to  reproduce  the  Kronecker  delta 
function,  all  the  frequencies  would  be  able  to  be  reproduced.  Unfortunately,  without  additional 
tools  this  is  not  easily  possible  without  large  number  of  particles,  due  to  the  stability  limit  of  the 
window  function  itself  requiring  a  minimum  of  two  nodes  in  its  support,  section  5.1. 

This  reconstruction  develops  confidence  in  the  method  and  helps  to  illustrate  the 
intricacies.  Most  notably,  it  will  show  the  difference  between  the  reproducing  kernel  method  and 
SPH  methods.  The  big  difference  is  the  effect  of  the  correction  function,  which  enables  an 
accurate  approximation  throughout  the  entire  domain  of  the  response  as  well  as  an  increased 
range  of  stable  operation,  especially  for  a  small  number  of  particles.  The  accuracy  of  the 
correction  function  can  be  analyzed  using  a  simple  reconstruction  of  a  known  function. 

It  is  readily  seen  in  figure  2  that  the  SPH  solution  is  not  even  able  to  reproduce  a  simple 
sinusoidal  wave  near  the  boundaries. 

In  order  to  relate  a  feel  for  the  performance  of  the  dilation  parameter,  figures  3  and  4 
depict  a  sinusoidal  reconstruction  with  and  without  the  correction  function.  Figure  3,  which  does 
not  contain  the  boundary  correction  function,  is  a  measure  of  the  SPH  method.  The  instability  of 
the  SPH  method  proves  that  the  GRK  correction  function  is  enhancing  stability  of  the  solution  as 
well  as  correcting  the  reconstruction  near  the  boundaries. 


Figure  3  {Left}  Sin(x)  reconstruction  using  the  GRK,  21  Nodes. 
Figure  4  (RightjSPH  reconstruction  for  sin(x)  using  21  nodes. 


II 


Figure  5  {Left}  Sin(x)  reconstruction  using  the  GRK,  201  Nodes. 

Figure  6  {Right}  SPH  reconstruction  for  sin(x)  using  201  nodes. 

It  is  important  to  note  that  j  =  -2  was  not  plotted  in  figure  4  since  the  instability  of  SPH 
method  would  have  required  are  much  larger  scale.  This  occurs  because  the  SPH  kernel  does  not 
satisfy  the  consistency  condition  for  small  number  of  nodes;  however,  for  large  numbers  of 
nodes  the  SPH  method  performs  adequately  (figure  6). 

Knowing  the  behavior  or  the  dilation  parameter,  namely  the  effect  of  the  j  parameter,  our 
discussion  is  turned  toward  a  brief  look  at  the  frequency  content  of  the  window  function.  By 
dilating  the  Gaussian  function  (figure  1),  it  is  easily  seen  how  the  function  changes  shape  with 
the  parameter  j.  It  is  this  dilation  that  enables  the  types  of  analyses  discussed  in  section  2.6. 

A  reconstruction  is  now  performed  on  a  sinusoidal  wave  that  has  been  augmented  to 
increase  its  frequency  as  well  as  dampen  its  amplitude. 


fix)  =  sin  ^^{x  +  .6  e  -0.1  jr  (3.0.1) 

The  results  in  figure  7  clearly  show  the  ability  of  the  window  function  to  capture  more  or  less  of 
the  frequency  content.  The  larger  the  values  for  the  refinement  parameter  the  more  quickly  the 
reconstruction  attempt  fails. 


Figure  7  The  reconstruction  of  a  sine  wave  with  increasing  frequency  and  damped  oscillations 

for  several  values  of  j  with  201  nodes. 
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4.0  Galerkin  GRK  Formulation 


The  reproducing  kernel  can  be  implemented  into  a  Galerkin  formulation  in  a  similar  way 
as  typical  finite  elements.  The  major  difference  in  construction  is  the  loop  occurs  over  nodes 
instead  of  elements,  but  the  formulation  is  almost  identical  beginning  from  Ae  weak  form  of  the 
momentum  equation.  The  following  variables  are  used:  t  is  the  traction,  b  is  the  body  force,  s  is 
the  stress  tensor,  and  u  denotes  acceleration. 


pSui'ui  dRx 


SuijGij  dRx 


JSuibi  dRx  +  I 
Rr  JdRt 


Suiti  dPx 


(4.0.1) 


This  paper  analyzes  a  one  dimensional  bar  with  a  step  force  applied  at  one  end  (figure  10).  The 
element  matrices  are  simple,  included  a  lumped  mass  matrix,  but  it  is  important  to  keep  in  mind 
that  the  matrices  have  a  variable  connectivity  depending  on  the  number  of  nodes  in  the  support  of 
the  window  function,  NEN. 


a,b  =  1. 


NEN 


(4.0.2) 


r'{x,  ay  =  =  j  I  (f  Na.x  0=1 .  NEN 


I 


(4.0.3) 


f\x,by  =  {<f^yj  =  ll  Nab^dR}<  a=l .  NEN 


■I 


I 


(4.0.4) 


r(x,  ty  =  ((/•^"l  =1  Nat^  dEx]  a  =  1 .  NEN 

\JdRxndRx 


(4.0.5) 


It  is  also  difficult  to  implement  the  boundary  conditions  since  several  shape  functions  can  be 
present  at  the  node  dictating  the  necessary  response.  The  natural  boundary  conditions  are  simply 
entered  into  the  force  vector;  however,  the  essential  boundary  conditions  are  more  difficult  to 
satisfy.  In  this  analysis  the  essential  boundary  condition  is  satisfied  by  substituting  a  sum  on  all 
the  shape  functions  in  place  of  the  row  in  the  stiffness  matrix  that  corresponds  to  the  constrained 
nodes.  The  shape  function  and  its  derivative  are  derived  with  the  correction  function  in  the  next 
two  sections. 


4.1  Properties  of  the  Shape  Functions 

Using  the  final  form  of  the  continuous  reproducing  kernel  (found  in  section  2.4),  it  is 
quite  easy  to  obtain  the  familiar  form  of  a  ualerkin  approximation.  If  the  equation  is  now 
numerically  integrated  and  terms  are  grouped  to  follow  the  standard  technique,  the  shape 


NP 


uHx)  =  X  Njix)  Uj 

(4.1.1) 

7=1 

NP 

u(x)=%  [c,  +  c4^^)] 

(4.1.2) 

7=1 

Nj(x)  =[  Ci(;c)  +  Ciix)  (^)]  ri  0{^^)aMj 

(4.1.3) 

The  characteristics  of  this  new  shape  function  must  be  carefully  analyzed  to  avoid 
erroneous  results.  The  most  apparent  difference  is  that  the  shape  function  does  not  meet  the 
Kronecker  delta  identity  since  each  node  is  influenced  by  several  shape  functions.  However,  Ae 
shape  function  will  meet  the  consistency  condition.  This  can  be  proven  by  the  following  equation 
(4.1.4)  where  the  reconstruction  equation  can  be  substituted  to  reproduce  itself.  It  is  also  pointed 
out  that  the  integration  method  used  to  calculate  C  and  C  -1  must  be  similar  in  order  to  cancel 
and  meet  this  condition.  It  is  also  proposed  by  Liu[93]  that  by  using  the  Trapezoidal  Rule  to 
integrate  the  moments  defined  within  C,  the  stability  of  the  kernel  is  increased. 


P(jc)=|  pPix)c\x)P'^is)  0i^-^)Pis)ds 
JRx 


=  P(Ar)rnx)j  p  P'^is)Pis)^  0i^)ds 
Jrx 

=  P(x)  (4.1.4) 

The  shape  functions  will  meet  the  following  isoparametric  shape  function  properties. 


Nj(x)  =  1 


NP 

and^  Nj{x)xj  =  x 


(4.1.5) 


4.2  Shape  Function  Derivatives 

Unlike  the  Moving  Least  Square  Methods,  Diffuse  Element  Methods,  and  Element  Free 
Galerkin  Methods,  the  derivative  of  the  GRK  shape  function  can  simply  be  obtained  by 
differentiation.  It  is  necessary  to  consider  the  correction  terms  as  well  as  the  window  function 
itself  to  obtain  the  derivative  of  the  shape  function. 


NjM  =  [c\{x)  +  C  2(x)  (^)  +  ^  j  r-i AMj 
+  (  C  i(x)  +  C  2(x)  (^) )  ’(^j  AMj 


(4.2.1) 


(4.2.4) 


C,  = 


mj  mi  (  mp  mu  +  mp  -  2m\  nti  ) 

momii-mi^  (momn-m^)^ 


Now  it  is  necessary  to  calculate  the  derivatives  of  the  moments  which  can  be  combined  with  the 
derivative  of  the  Gaussian  window  function  to  reveal  the  final  derivative  in  equation  (4.2.8). 


m'o(x)=l  /Y  ds  (4.2.5) 

m'  i(x)  =  I  ^^[2  -  1]  ds  (4.2.6) 

Jbm’' 

m'  ii(a:)  =  f  [(^F  -  1]  e-Y^Y  ds  (4.2.7) 

Jb(x) 

NjAx)  =  j  C  ;(x)  +  C  2ix)  (^)  +  ^  j  r  -1  AMj 

^  (4.2.8) 

-^(Ci(x)  +  C2(x)(^))(^)r-i  0(^)  AMj 


5.0  Stability  Analysis 

It  is  interesting  to  note  that  the  stability  of  this  type  of  analysis  is  two-fold.  First  there  is 
the  stability  of  the  Reproducing  Kernel  itself,  and  then  there  is  the  stability  of  the  time 
integration  method.  In  this  analysis  an  Explicit  Newmark  Beta  Predictor/Corrector  Type 
algorithm  is  implemented  (Hughes  and  Liu  [78]). 

5.1  Reproducing  Kernel  Stability 

The  stability  of  the  kernel  is  mainly  a  function  of  the  number  of  nodes  encompassed  by 
the  Gaussian  window  function.  (Theoretically  the  number  of  nodes  covered  is  the  number  of 
nodes  in  the  analysis,  since  the  function  has  an  infinite  domain.)  The  number  of  significant  nodes 
in  a  shape  function  is  controlled  by  the  refinement  parameter,  r.  From  the  aforementioned 
definitions  for  the  shape  of  the  window  function  the  following  equation  can  be  derived  to 
estimate  the  number  of  nodes  under  any  given  window  function.  If  the  radius  of  a  given  window 
function,  Dxc,  is  defined  as  the  distance  from  the  center  to  the  edge  of  the  windows  significant 
support . 

Then  a  ratio  can  be  defined  to  relate  the  height  of  the  Gaussian  function  at  the  peak  to  the 
small  value  where  it  is  safely  approximated  to  be  zero. 
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(5.1.2) 


(• 


Using  these  definitions  it  is  possible  to  derive  the  kernel  stability  stated  in  equation  (5.1.3)  where 
the  variables  are  as  follows. 

n  =  the  number  of  nodes  covered  by  the  Gaussian  function 
j  =  the  parameter  controlling  the  dilation 

R  =  the  ratio  between  the  value  at  a  node  in  the  tails  of  the  function  to  the 
significance  of  the  peak  of  the  Gaussian  function. 

+  (5.1.3) 

Theoretically  the  number  of  nodes  should  be  at  least  two  in  order  to  maintain  the  variable 
connectivity  window  arrays.  This  is  verified  analytically  through  solutions  that  were  run  with 
values  as  low  as  -2.5  for  j  when  the  entire  domain  is  used  as  the  support  of  the  Gaussian  function, 
and  as  low  as  -2.2  when  the  support  of  the  Gaussian  function  is  limited.  Solutions  for  the  number 
of  nodes  in  the  support  of  a  window  function  for  various  height  ratios  are  shown  in  the  following 


j 

n  (2s) 

‘n  (3s) 

-2.2 

1.69 

2.04 

-2 

1.80 

-1 

2.60 

3.39 

0 

4.19 

5.79 

1 

7.38 

10.57 

2 

13.77 

20.15 

3 

26.53 

39.30 

NOTE;  2s  and  3s  corresponds  to  R  =  1.83e-2  and  R  =  1.23e-4  respectively. 

5.2  Critical  Time  Step 

In  order  to  perform  analyses  on  the  structural  dynamic  class  of  problems,  it  is  very 
important  to  understand  the  relationship  between  this  new  shape  function  and  the  critical  time 
step.  This  is  determined  here  by  performing  the  standard  eigenvalue  analysis,  solving  the 
determinant  of  [K  - 1 M].  This  calculation  only  needs  to  be  performed  for  the  element  case  since 
it  is  known  that  the  maximum  element  frequency  determines  the  upper  bound  for  the  critical  time 
step.  The  critical  time  step  was  found  to  depend  on  several  parameters  including  the  dilation 
parameter  and  the  boundary  correction  function. 

Using  symbolic  manipulation  this  determinant  was  solved  for  GRK  shape  functions 
containing  several  different  numbers  of  nodes.  By  evaluating  the  determinant  it  was  found  that 
there  is  only  one  non-zero  eigenvalue,  regardless  of  the  number  of  nodes  in  the  support  of  the 
shape  function.  The  result  was  able  to  be  simplified  for  equally  spaced  nodes  to  equation  (5.2.3) 
by  defining  the  following  terms. 


£)(x,j)=  Ci(,x)  + C2(x)(^ 

(5.2.1) 

D'  ix,s)  =  C  'i(x)  +  C  2(x)  (^)  -h  ^ 

(5.2.2) 
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=  [2^] dx^  [a  <P.,-  +  D'i  (5.2.3) 

where 

dx  is  the  length  between  nodes  (density  =  1) 

Dxi  is  the  integration  weight  using  Trapezoidal  Rule 
NP  is  the  number  of  nodes  in  the  shape  function’s  support 
Finally,  the  critical  time  step  is  shown  for  a  central  difference  time  integration  scheme  in 
equation  (5.2.4).  By  substituting  the  maximum  eigenvalue  corresponding  to  the  maximum 
frequency,  the  critical  time  step  can  be  calculated. 

At<  (5.2.4) 


-1/2 

(5.2.5) 


Results  for  the  critical  time  step  were  calculated  using  equation  (5.2.5).  The  results  are 
shown  in  the  following  table  for  2 1  nodes.  This  can  be  compared  to  the  standard  finite  element 

At  <1/  c  which  is  2.455  x  lO"^  s. 


At<  ^^^dx 


NP 


Y[Di<!>'.i+D'i(!>,iYAxi 


j 

At  (Analytic) 

At  (Numeric) 

-1 

0 

1 

4.912  X  10-^ 

5.651  X  10-5 

The  discretized  form  of  the  critical  time  step  was  perceived  to  have  an  simple  translation 
to  the  continuous  form.  The  discretized  form  contains  the  derivative  of  the  reproducing  kernel 
squared.  Then  the  largest  eigenvalue  and  thus  the  critical  time  step  is  as  follows: 


A  =|^j  r2  f  [D(x,y)  +  D\x,y)  dy 

A  =  1^-^]  I  dy 

'  ^  '  JB(f) 

where  k(x,y)  is  the  reproducing  kernel  defined  by  equation  (5.2.8). 
/:(x,y)  =  D(x,y)<^(^) 


(5.2.6) 


(5.2.7) 


(5.2.8) 


6.0  Numerical  Experiments 
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steel)  E  =  3  X  lO^  psi,  r  =  7.24  x  lO'^  slugs/in^,  and  A  =  1  in^.  For  the  elastic  -  plastic  problen, 
yield  stress  ay  =  30000  psi,  Ep  =  E/4  and  Fq  =  5000  lb  are  employed. 


- h-  Fo  =  1000  lb 

Length  =  100" - - 


Figure  8  1-D  linear  elastic  rod  with  step  input. 


The  solution  for  the  linear  elastic  deforming  rod  was  obtained  with  the  GRK  and  a 
Explicit  Newmark-Beta  predictor/corrector  algorithm.  There  results  coincided  accurately  with 
the  closed  form  solution. 


Time  (s) 


Figure  9  Displacement  for  node  1 1  (center  node)  showing  a  range  for  the  dilation  parameter 

and  the  exact  solution. 

In  order  to  visualize  the  effect  of  the  dilation  parameter  in  this  application,  the  shape 
function  are  plotted  in  figure  10  along  with  the  nodes  which  depicts  the  number  of  significant 
nodes  providing  the  support  of  each  shape  function  in  each  dilation  (Technically  every  node  is  in 
the  support  of  every  Gaussian  shape  function.)  This  is  shown  directly  above  the  plot  for  velocity 
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Velocity  (in/s) 


versus  time  which  incorporated  the  shape  functions  in  figure  10  to  obtain  this  response  (figure 

11). 


Figure  10  Shape  functions  used  to  calculate  the  wave  propagation,  with  nodal  support. 


Figure  11  Velocity  for  node  1 1  (center  node)  showing  a  range  for  the  dilation  parameter  and 

the  exact  solution. 


For  completeness  the  stress  is  also  plotted,  figure  14. 


Figure  12  The  axial  stress  is  plotted  for  the  middle  node  versus  time, 


Figure  13  Axial  stress  at  the  tenth  time  step  as  a  function  of  position, 


Time  (s) 

Figure  14  GRK  run  with  j  =  0  and  a  time  step  of  3.32  x  10’^  which  is  1.35  times  the  FEM 

critical  step. 


41  nodes,  de«  =  0.00001 


time 


Figure  15  Elastic-plastic  responses  at  the  midpoint  of  the  rod. 

Although  the  solution  in  figure  14  is  not  as  accurate  as  the  finite  element  solution,  it  is 
interesting  to  note  that  the  problem  was  run  at  1.35  times  the  FEM  critical  time  step. 
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7.0  Conclusions 

The  theory  for  a  Gaussian  Reproducing  Kernel  is  presented  here  with  numerical 
experiments  performed  to  confirm  the  derived  equations.  The  GRK  method  is  proven  here  to 
possess  the  ability  to  solve  a  dynamic  problem.  Results  are  also  presented  to  verify  the 
supposition  that  the  correction  function  can  provide  both  boundary  correction  and  reproducing 
kernel  stability.  Furthermore  the  implementation  of  the  flexible  window  function  as  frequency 
control  has  been  initiated. 

The  results  from  the  numerical  experiments  not  only  verified  the  theory  presented,  but  it 
produced  several  encouraging  results.  Among  the  most  important  is  the  ability  of  the  Gaussian 
Reproducing  Kernel  to  perform  at  time  steps  larger  than  the  critical  time  step  for  standard  finite 
elements.  It  is  also  important  that  results  for  the  correction  function  proved  that  it  provided 
boundary  correction  as  well  as  enhancing  the  stability  of  the  reconstruction  equation.  The 
increased  stability  of  the  reconstruction  equation  enables  the  ability  to  use  fewer  particles  in 
analyses.  The  correction  function  increased  the  stable  range  of  the  dilation  parameter  over  the 
SPH  method. 

The  advantageous  characteristics  of  the  GRK  were  not  able  to  produce  superior  results 
over  the  standard  finite  element  method  due  to  the  finite  element’s  inherent  ability  to  solve  this 
elastic-dynamic  problem  exactly  when  run  at  the  critical  time  step.  Nonetheless,  the  theory  has 
been  derived  and  it  is  logical  that  characteristics  such  as  continuous  derivatives,  flexible  window 
filtering,  and  the  kernel’s  increased  stability,  and  increased  critical  time  step  ability  are  valuable 
tools  that  can  be  exploited  in  further  experiments. 
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